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SUMMARY 


The influence of changes in precursor region flow properties 
(resulting from the absorption of radiation from the shock layer) 
on the entire shock-layer flow phenomena is investigated. The 
axially symmetric case is considered for both the preheating zone 
(precursor region) and shock layer. The flow in the shock layer 
is considered to be viscous with chemical equilibrium but radiative 
nonequilibrium. Realistic thermophysical and spectral models are 
employed, and results are obtained by implicit finite difference 
and iterative procedures. The results indicate that precursor 
heating increases the radiative heating of the body by a maximum 
of 7.5 percent for 116-km entry conditions. 
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1 . 


INTRODUCTION 


Radiation plays a very important role in the analyses of flow 
phenomena around an entry body at high-speed entry conditions. In 
many instances, the radiative energy transferred to the body exceeds 
the convective and aerodynamic heat transfers (refs. 1-10). 

Radiative energy transfer from the shock layer of a blunt body 
into the free stream reduces the total enthalpy of the shock layer 
while increasing the enthalpy of the free-stream gases. Because 
of this increase in enthalpy, the entire flow field ahead of the 
shock layer and around the body is influenced significantly. The 
precursor flow region is considered to be the region ahead of a 
shock wave in which the flow field parameters have been changed 
from free-stream conditions due to absorption of radiation from 
the incandescent shock layer. Most of the radiative energy trans- 
ferred from the shock layer into the cold region ahead of the 
shock is lost to infinity unless it is equal to or greater than 
the energy required for dissociation of the cold gas. When the 
photon energy is greater than the dissociation energy, it is 
strongly absorbed by the cold gas in the ultraviolet continuum 
range. The absorbed energy dissociates and ionizes the gas, and 
this results in a change of flow properties in the precursor 
region. In particular, the temperature and pressure of the gas 
is increased while velocity is decreased. The change in flow 
properties of the precursor region, in turn, influences the flow 
characteristics within the shock layer itself. The problem, 
therefore, becomes a coupled one, and iterative methods are 
required for its solution. 

Only a limited number of analyses on radiation-induced pre- 
cursor flow are available in the literature. Works available 
until 1968 are discussed, in detail, by Smith (ref. 11) . By 
employing the linearized theory of aerodynamics. Smith investi- 
gated the flow in the precursor region of a reentry body in the 
Earth's atmosphere. The cases of plane, spherical, and cylin- 
drical point sources were considered, and solutions were obtained 
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for a range of altitudes and free-stream conditions. It was 
found that for velocities exceeding 18 km/s, precursor flow 
effects are greatest at altitudes between 30 and 46 km. It 
was further concluded that preheating of air may cause an 
order of magnitude increase in the static pressure and tempera- 
ture ahead of the shock wave for velocities exceeding 15 km/s . 

In a preliminary study, Tiwari and Szema (ref. 12) investigated 
the change in flow properties ahead of the bow shock of a 
Jovian entry body resulting from the absorption of radiation 
from the shock layer. The analysis was done by employing the 
small perturbation technique of classical aerodynamics as well 
as the thin layer approximation for the precursor region. By 
employing appropriate thermodynamic and spectral data for the 
hydrogen/helixom atmosphere, variations in precursor region flow 
quantities (velocity, pressure, density, temperature, and enthalpy) 
were calculated by the two entirely different methods. For Jovian 
entry conditions, one-dimensional results obtained by the two 
methods were found to be in good agreement for the range of 
parameters considered. It was found that preheating of the gas 
significantly increases the static pressure and temperature ahead 
of the shock for entry velocities exceeding 35 km/s. It was con- 
cluded that for certain combinations of entry speeds and altitudes 
of entry, the precursor effects cannot be ignored while analyzing 
flows around Jovian entry bodies. Specifically, it was seen that 
at an altitude of 95 km, the precursor effects are important for 
entry velocities greater than 35 km/s. 

In the analyses of most shock-layer flow phenomena, the con- 
tribution of radiation-induced precursor effects usually is neg- 
lected. However, a limited niomber of analyses which include this 
effect are available in the literature. Lasher and Wilson (refs. 
13, 14) investigated the level of precursor absorption and its 
resultant effect on surface radiation heating for the Earth's 
entry conditions. They concluded that, for velocities less than 
18 km/s, precursor heating effects are relatively unimportant in 
determining the radiative flux reaching the surface. At velocities 
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greater than 18 km/s, the amount of energy loss from the shock 
layer and resultant precursor-heating correction was found to 
be significantly large. Liu (refs. 15, 16) also investigated 
the influence of upstream absorption by cold air on the stag- 
nation region, shock layer radiation. The thin layer approxi- 
mation was applied to both the shock layer and the preheating 
zone (the precursor region) . The problem was formulated for 
the inviscid flow over smooth blunt bodies, but the detailed 
calculations were carried out only for the stagnation region. 

The general results were compared with results of two approxi- 
mate formulations. The first approximate formulation neglects 
the upstream influence, and the second one essentially uses the 
iterative procedure described by Lasher and Wilson (refs. 13, 

14) . The results are compared for different values of a 
radiation/convection parameter. A few other works, related to 
the effects of upstream absorption by air on the shock-layer 
radiation, are discussed by Liu (refs. 15, 16) . Some works on 
precursor ionization for air as well as hydrogen/helium 
atmosphere are presented in references 17 through 21. 

The first purpose of this study is to investigate the flow 
properties in the entire precursor region ahead of a Jovian 
entry body. To accomplish this, the precursor region is assumed 
to be thin, and the flow in this region is considered to be 
inviscid. In this respect, therefore, the proposed study may 
appear as an extension of the analysis presented in reference 12. 
Next, it is proposed to investigate the influence of changes in 
the precursor region flow properties on the entire shock-layer 
flow phenomena. The flow in the shock layer is considered to be 
viscous and in chemical equilibrium. The solutions of the 
governing viscous shock-layer equations are obtained by employing 
the numerical procedures outlined in reference 9 and 22. It 
should be emphasized again that the flow phenomena in the shock 
and precursor regions are coupled, and iterative procedures are 
needed for finding solutions. In this respect the proposed 
investigation differs from the analysis presented in reference 12. 
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The basic formulation of the problem is presented in section 
2. Radiation absorption models and radiative flux equations are 
given in section 3. The data required for finding the solutions 
of the governing equations are given in section 4 , and the 

solution methods are presented in section 5. The results are 
obtained for the Jupiter ' s entry conditions and for an entry 
body which is a 45® hyperboloid. Precursor as well as shock- 
layer results are discussed in section 6. 

2. BASIC FORMULATION 

The physical model and coordinate system for a Jovian entry 
body are shown in figure 2.1. The entire flow field ahead of the 
body can be divided essentially into three regions: the free 

stream, the precursor region, and the shock layer. The flow 
properties are considered to be uniform at large distances from 
the body. In this section, basic governing equations and the 
boundary conditions are presented for the precursor as well as 
shock-layer regions. 


2.1. Precursor Region 

In this region, the flow is considered to be steady and 
inviscid. With reference to the coordinate system shown in figure 
2.1, the governing equations for an axisyiranetric flow can be 
written as (refs. 12, 23) 

Mass continuity: 

(9/3s) (pur) + (3/3n) (pvXr) = 0 (2.1) 

Momentum: 

p[u(9u/9s) + Xv(9u/9n) - Kuv] + (9p/3s) = 0 (2.2) 

p[u(9v/3s) + v(9v/3n) + Ku^] + X(3p/3n) = 0 (2.3) 
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Energy : 


p [ (u/X) (9H/8S) + v(9H/3n)] + (Xr) ^ [ ( 9/ 9n) (Xrq^^) ] = 0 (2.4) 

Species continuity : 

p[(u/5i) (ac^8s) + vOC^an)] - = O (2.5) 

where K = K(s) = l/R^, X = 1 + Kn. It should be noted that, 
according to the notations used in figure 2.1, all quantities 
appearing in the above equations should have a prime superscript 
(i.e., u', v', p', H', etc.), and all physical coordinates 

should have a superscript * (i.e., s*, n*, r*, etc.). 

However, for the sake of clarity, these notations have been 
omitted from the equations. 

If the precursor region is assumed to be thin, then one can 

make the approximation that (n/R ) << 1, (9/9s) , (9/9n) , and R 

s s 

is not a function of n. In this case, X = 1, and equations 
(2.1) through (2.5) reduce to 


(9/9n) (pv) 

= 0 

(2.6) 

pv(9u/9n) = 

0 

(2.7) 

pv(9v/9n) + 

(9p/9n) == 0 

(2.8) 

pv(9H/9n) + 

(9qj^/9n) = 0 

(2.9) 

pv(9c^/9n) - 

K =0 
a 

(2.10) 


The boundary conditions for this region are the free-stream 
conditions and the conditions at the outer edge of the shock. 

For the coupled precursor/shock- layer flow phenomena, the boundary 
conditions at the outer edge of the shock are obtained through 
iterative procedures. 
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2.2. Shock-Layer Region 


In this region, the conditions for which the present analysis 
is carried out are that the flow is axisyinmetric , steady, laminar, 
and compressible. It is further assumed that the gas is in the 
local thermodynamic and chemical equilibrium, and that the tangent 
slab approximation is valid for radiative transport. For this 
region, the viscous shock-layer equations presented in references 
9 and 22 are a set of equations that are valid uniformly through- 
out the shock layer. The methods of obtaining these equations are 
discussed in detail in the those references. First the conservation 
equations are written for both the inviscid and the boundary- layer 
regions in the body-oriented coordinate system. Then these equa- 
tions are nondimensionalized in each of the two flow regions with 
variables which are of order one. Terms in the resulting sets of 
equations are retained up to second order in the inverse square 
root of Reynolds number. Upon combining these two sets of 
equations so that terms up to second order in both regions are 
retained, a set of equations uniformly valid to second order in 
the entire shock layer is obtained. The nondimensional form of 
the viscous shock-layer equations that are applicable in the 
present case can be written as 

Continuity : 

(3/3x) |^(r + y cos 0)pu + (3/3y) (1 + y<) (r + y 

X cos 0)pvj =0 (2.11) 

x-momentum: 

/ u 3u 3^ uvK \ 1 3p 

^\1 + y< 3x ^ 9y 1 + yK / 1 + yK 3x 



( 2 . 12 ) 
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y-momentum ; 


Energy ; 


( ^ 

9v 

Vi + yK 

!)x 

( ^ 

9H 

M + yK 

9x 

: e4A[jL . 
\0yLpr 

]JKU^ 

.1 . 


(2.13) 


N 


N. 


COS 


+ Ji- (Pr - Du |a 


Pr 


_0 \[~ u 8H __ JJ_ 


N 


1 + yKj V 1 + yK r + y cos 0 9y Pr 

- g, Vi - If - T^Jl - [ 


9C. 


3_5r 

9y 


+ q 


R 


( 1 

\1 + 


COS 0 


yK 


)] 


(2.14) 


where H = h + u^/2. 

The terms used to nondimensionalize the above equations are defined 
as 


X = x*/R* 
n 

y = y*/R* 

n 

r = r*/R* 

' n 

■3r = 

h = h*/V*^ 

J. = J*R*/u* ^ 
1 in ^ref 


V = v*/V* 

OO 

p = p*/p^; 


U = P*/Pjef 


= c*/c* 

P P Poo 


^ref = 


K = K*/R* 
n 




(2.15) 
(cont ' d. ) 
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P = pV(p*v*") 


Pr = C*y*/K* 


T = T*C /V*^ 

pco' 00 


Le. . = p*C*D* /K* 
ID ^ P ID 


u = u*/V^ 


L. . = o*C*D* ./K* 
ID ^ P 13 


(2.15) 
(concl ' d . ) 


In equation (2.14) , J. represents the mass flux relative 
to the mass average velocity and is given by the expression 
(refs. 24, 25) 


J. 

1 


= - (P/Pr) 


_ m 

E (l!/T) OT/3y) 

k“ 1 


(2.16a) 


where 


^iK 


Le^, i = K 
i K 


NI 




NI 




j?^i 


NI 

Ab^j^ = Le^ - |(M^/M)Le^j^ + [1-(M^/I^)] .L^be^jCj 


The last term in equation (16a) represents the contribution of 
thermal diffusion. The quantity be^^ represents the multi- 
component Lewis number, and represents the binary Lewis 

Semenov numbers; both are defined in equation (2.15). If thermal 
diffusion can be neglected and L^^ can be taken as constant 
for all species, then equation (2.16a) reduces to 
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J . 

1 


(y/Pr)L. . OC./3y) 

-i*J ^ 


(2.16b) 


In the present study, use is made of equation (2.16b) , and the 
value for is taken to be 1.1 (ref. 26). 

The expression for the equation of state for a hydrogen/ 
helium mixture is given by Zobi et al. (ref. 27) as 

T* = [ (p*/1013250)V (p*/0. 001292)^] (2.17a) 


H* = C„ [ (p*/1013250)^/(p*/0. 001292)^1 (RT /M) (2,17b) 

ti o 

where 


0.65206 - 

0.04407 

in 

0.67389 - 

0.04637 

in(Xn^ 

0.95252 - 

0.1447 

(Xh, 

0.97556 - 

0.16149 

«n(Xn, 


U. = sin 0[1 + 0,7476(1-X„ )] 

T. oo H2 

CTU = —545.37 +• 61.608 — 22459 + 0.039922 

- 0.00035148 + 0.0000012361 U® 

CHU - 5.6611 - 0.52661 U. + 0.020376 U? - 0.00037861 U? 

t t t 

+ 0.0000034265 - 0.000000012206 U® 

^ t 

= CTU + 61.2(1 - Xjj^) 


C„ = CHU - 0.3167(1 - X„ ) 
H Ha 
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and X„ represents the mole fraction of H 2 . 

112 

The set of governing equations presented above has a 
hyperbolic/parabolic nature. The hyperbolic nature enters 
through the normal momentum equation. If the shock layer is 
assumed to be thin, then the normal momentum equation can be 
expressed as 


Pu^k/(1 + y<) = Op/3y) 


(2.18) 


If equation (2.13) is replaced with equation (2.18), then the 
resulting set of equations is parabolic. These equations can, 
therefore, be solved by using numerical procedures similar to 
those used in solving boundary-layer problems (refs. 9, 22). 

In order to solve the above set of governing equations, it 
is essential to specify appropriate boundary conditions at the 
body surface and at the shock. At the body surface, the no-slip 
boundary conditions are used in this study. Thus, the following 
conditions are imposed at the body surface: 

u = V = 0 (2.19) 

T = Tj^ = constant (2.20) 

The conditions in front of the shock are obtained from the 
solution of the precursor region flow field. The conditions 
immediately behind the shock are obtained by using the Rankin- 
Hugoniot relations. The nondimens ional form of the shock 
relations are expressed as 

Mass continuity: 

P„- v' = " sin a (2.21) 

w O 
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Momentum : 



cos a 


( 2 . 22 ) 


p = p + + sin^ a[l - (1/p ..)] 

o S ^ 


(2.23) 


Energy: 


hg_ = hg+ + [sinMa/2)][l - (1/Pg-)1 (2.24) 

Note that in the above equations u' and v' represent the 
velocity components in the shock-oriented coordinate system. 

The relations for u^ and v^ in the body-oriented coordinate 
system can be obtained from 

u = u' sin(a + 3) + v' cos(a + 3) (2.25) 

s s s 

V = - u' cos (a +3) + v' sin(a + 3) (2.26) 

s s s 

The governing equations and the boundary conditions pre- 
sented in the above two subsections essentially describe the 
flow field in the precursor/shock-layer regions. In the next 
two sections appropriate relations for the radiative flux, 
thermodynamic and transport properties, and equilibrium composi- 
tion of the gas will be presented. 

3, RADIATION MODELS 

An appropriate expression for the radiative flux q^^ is 
needed for the solution of the energy equation presented in the 
previous section. This requires a suitable transport model and 
a meaningful spectral model for variation of the absorption 
coefficient of the gas. In this section, appropriate expressions 
for the spectral and total radiative flux are given and informa- 
tion on the spectral absorption by the hydrogen/helium gas is 
presented. 
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3.1. Radiative Flux Equations 

The equations for radiative transport, in general, are 
integral equations which involve integration over both frequency 
spectrum and physical coordinates. In many physically realistic 
problems, the complexity of the three-dimensional radiative 
transfer can be reduced by introduction of the "tangent slab 
approximation." This approximation treats the gas layer as a 
one-dimensional slab in calculation of the radiative transport. 
Radiation in directions other than normal to either the body or 
shock is neglected in comparison. Discussions on the validity 
of this approximation for planetary entry conditions are given 
in references 28 to 32. 

As mentioned earlier, the tangent slab approximation for 
radiative transfer is used in this study. It should be pointed 
out here that the tangent slab approximation is used only for 
the radiative transport and not for other flow variables. For 
a nonscattering medium and diffuse nonreflecting bounding surfaces, 
a one-dimensional expression for the spectral radiative flux is 
given by (refs. 33, 34): 

= 2,rje^[B^(0)E,(T^) - - x^)l 

o 

where 

J (y* )dy ' 
o 

En(t) = J' exp(- t/y)y^ ^ dy 
o 

= (hv^/c^) [exp(hv/KT) - 1] 
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The quantities (0) and represent the radiosities 

of the body surface and shock respectively. 

The expression of total radiative flux is given by 
o 

To obtain specific relations for the total radiative flux for 
the precursor and shock-layer regions, it is essential to know 
the spectral absorption characteristics of the absorbing-emitting 
species in these regions. 

In the precursor region, the radiative contribution from 
the free stream usually is neglected. For a diffuse, nonreflecting 
shock front, the expression for one-dimensional radiative flux 
for this region is obtained from equations (3.1) and (3.2) as 

00 

qj^(n) = 2f {q^ (0) E3 (K^n) 
o 

00 

+ TTK^ r (T)E 2 [k (n - n')]dn'}dv (3.3) 

V V V 

where q^(0) = e^'iTB^(Tg). In obtaining the above equation it was 
assumed that the absorption coefficient is independent of 

position. 

The information on the spectral absorption model for hydrogen/ 
helium species in the precursor region is given in reference 12 
and is briefly discussed in subsection 3.2. The model essentially 
consists of approximating the actual absorption of active species 
by three different step models. For this model, equation (3.3) 
can be expressed as (ref. 12) : 
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N 

q^in) = 2ir^ 
i=l 


V . 
21 


( 15 /TT®)q( 0 )E 3 (K^n) J' [vV(e^ - l)]dv 


^li 


+ K 


n 

i/ ' 


E 2 [k^ ( n - 


V 

n')l / 
Vli 


ii 


(T)d^dn' 


(3.4) 


where v = hv/KT^ and q(0) = £crT^. In writing the above equation 
it has been assumed that the shock front radiates in the pre- 
cursor zone as a gray body. 


In the shock layer, the radiative energy from the bow 
shock usually is neglected in comparison to the energy absorbed 
and emitted by the gas layer. This implies that the transparent 
shock front does not absorb but emits radiation. The expression 
for the net radiative flux in the shock layer, therefore, is given 
by 


= 2 / 


q^(0)E3(T^) +/ ^ B^(t)E2(x^ - t)dt 


ov 


-/ B (t)E2 (t - T )dt 


dv 


(3.5a) 


In this equation, the first two terms on the right represent 

the radiative energy transfer towards the bow shock while the 

third term represents the energy transfer towards the body. 

Upon denoting these contributions by q and q , equation 

R R 

(3.5a) can be written as 



A few spectral models for absorption by the hydrogen/ 
helium species in the shock layer have been proposed in the 
literature (refs. 35 to 37). For Jovian entry conditions. 
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the absorption by helium usually is neglected. The spectral 
absorption of hydrogen species was represented by a 58-step 
model by Sutton (ref. 36) and was approximated by a 30-step 
model by Tiwari and Subramanian (ref. 37) . The results of these 
step models are compared with the detailed model of Nicolet (ref. 
35) in reference 37. The 58-step model proposed by Sutton is 
employed in this study. The details of radiative absorption and 
computational procedure are given in references 36 and 37. The 
information on spectral absorption by this model is summarized 
in subsection 3.2. In essence, the step model replaces the 
frequency integration in equation (3.5) by a summation over 
58 different frequency intervals. In each interval, the 
absorption coefficient is taken to be independent of frequency. 
For this model, equation (3.5) can be expressed as 


N 

= 2 -, Y. 

j=l 


{ 

' 

{e B (T )E3 
1 V V w ^ 

O _ 




s 


(y ' ) dy ' 


dc 



a^(?)B^(UE2 




(y ' )dy ' 


d? 


(3.6) 


where y^ denotes the shock location and N represents the 
nximber of spectral intervals. In each of the jth intervals, 
the absorption coefficient is assumed constant while the Planck 
function is not. In accordance with equation (3.5b), equation 
(3.6) can be expressed in terms of q!^ and q", and for a 
gray body one finds 


qR(y) 


N 

(4TTh/c2) ^ 

j=l 


eF (v . rT^^).Es 
J W 




V3 


(y')dy 


I 


+ (KT/h) 
o 


F(Vj,T)a^j (^)E2 



y 


V3 


(y ' ) dy 



(3.7a) 
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q;(y) 


N 

(4TT.h/c^) X) 

j=l 


I f (KT/h) 


F{Vj ,T)a^jE2 


/ 


“vj 


X 




(3.7b) 


where 


F(v.,T^) = J {v [exp (hv/KT^) - l]}dv 


''j' 


V . 


F(v.,T) = f ^ {v®/[exp(v) - l]}dv, V = 

^ T T 


hv/KT 


''j* 


From the knowledge of the temperature distribution normal 
to the body, equations (3.7) can be solved by numerical integration 
over frequency and space. The final temperature profile is 
obtained through an iterative procedure. Use of equations (3.7) 
is made in obtaining the radiative flux towards the body and shock 
as well as the net radiative flux. 

For evaluation of the radiative flux, usually it is essential 
to express the exponential integrals E^(t) in simpler approximate 
forms. Quite often, these integrals are approximated by 
appropriate exponential functions (refs. 33, 34). In this study 
it was established that better results are obtained if the 
exponential integrals are expressed in series form for small and 
large arguments. The series expansion of the exponential integral 
of first order is given as 


For t < 1: 


El (t) 


0.5772 


iln t + t 


t" 

2(2) I 


t 

3(3) 


+ . . . 


For t > 1 : 

El (t) 


exp (-t) 


ao + ait + aat^ + ast^ + t** 
t(bo + bit + hzt^ + bjt^ + t^) 


(3.8a) 


(3.8b) 
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where 


ao = 0.26777343 
ai = 8.63476089 
a2 = 18.059016973 
as = 8.57322874 


bo = 3.958469228 
bi = 21.09965309 
b2 = 25.63295614 
bs = 9.5733223454 


Relations for exponential integrals of higher order are obtained 
by employing the recursion relations given in reference (34). 

3.2. Radiation Absorption Model 

Appropriate spectral models for gaseous absorption are needed 
for solutions of the radiative flux equations. Information on 
spectral absorption by the precursor and shock-layer species is 
presented in this subsection. 

3.2.1, Spectral absorption model for precursor region : - In 
the precursor region, the photoionization absorption coefficient 
is a continuous nonzero function of photon energy (because of 
bound-free transition) for all values of photon energy that exceed 
the ionization potential of the atom. Similar remarks apply to 
the photodissociation and radiative recombination. A critical 
review of ultraviolet photoabsorption cross sections for molecules 
of astrophysical and aeronomic interest, available in the litera- 
ture up to 1971, is given by Hudson (ref. 38). Specific infor- 
mation on photoionization and absorption coefficient of molecular 
hydrogen is available in references 38 to 44. 

Photoionization and absorption cross sections of H 2 , as 
obtained from references 38 to 44, are plotted in figure 3.1. 

From this figure it is evident that the ionization continuum 

O 

starts at about 804 A and continues towards lower wavelengths. 
Between the wavelengths of 600 A and 804 A, the absorption cross 
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section for the ionization continmam is included in the total 
absorption (i.e., absorption due to ionization as well as 

O 

dissociation) . For wavelengths below 600 A, however, the 
ionization continuum absorption is equal to the total absorption. 
The total absorption cross section for the continuum range 

O 

below 804 A can be closely approximated by the two rectangles 
(I and II) shown in the figure with broken lines. The ratio of 
the ionization cross section to the total absorption cross- 
section (i.e., the value of Y^) is taken to be unity for 
rectangle I and 0.875 for rectangle II. For wavelengths greater 

O 

than 804 A (where hv is below ionization energy) , the value 
is taken to be zero. Little information is available in 
the literature on the absorption cross section for dissociation 
of H 2 molecules. There is strong evidence, however, that 

O 

photodissociation starts at about 2600 A and continues towards 

O 

lower wavelengths to about 750 A (refs. 39 to 41). There are 
also a few diffuse bands in this spectral range (refs. 39, 41). 
Thus, it becomes difficult to evaluate the absorption cross 
section in this spectral range. For this study, the absorption 

o o 

cross section in the spectral range between 804 A and 2600 A was 
approximated by rectangle III. The specific values of a(v) 
for the three rectangles are found to be a^(v) = 4.1 E-18, 
cTjj.(v) = 8.2 E-18, and cfj^j(v) = 2.1 E-18. The value of Y^^ is 
taken to be zero for rectangle I and 0.125 for rectangle II. The 
numerical procedure for employing this model in the radiative 
flux equations is discussed in detail in reference 12 and is 
summarized in section 5. 

3.2.2. Spectral absorption model for shock layer : - As 
mentioned earlier, the 58-step model proposed by Sutton (ref. 36) 
for spectral absorption by the hydrogen species in the shock 
layer is employed in this study. For atomic hydrogen, all three 
transitions (bound-bound, bound-free, and free-free) are 
considered. The total absorption of the jth step is a simrmation 
of the average absorption for the ith transitions in the jth 
step, i.e. , 
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(3.9b) 


K 


ij 


V .+Av . 

(1/Av.) ^ K^dv 


= f(T,N^,V) 


(3.9c) 


where represents the number density in cm ® . 

For the free-free transition, the absorption coefficient 
is calculated by 

Kff = (2.61E - 35)N^Njj+/(v*T^/^) (3.10) 


The absorption coefficient for bound-free transitions is 
calculated by employing two separate relations as 

(3.11a) 
(3.11b) 

Cl = (-157780/T) [1 - d/nj)3' 

C2 = (-157780/T) (1 - 6/13.6) 

Ca = [ (157780/T) (1/25 - 6/13.6)] - 1 

6 = (1.79E - 5) (N^/^)/(T^'^^) 

In the above equations, n^ represents the principal quantum 
numbers, 6 is the reduction in ionization potential in eV, 


H 

bf 


= (1.99E - 14) (Njj/v’) ^ (1/n’ ) exp(Ci), 


n^=l 




= (6.31E - 


20) (TN^/v^) exp(C 2 ) exp(C 3 ), 


5 < n„ < n. 

— Z — Jl,max 


where 



and the values 157780 and 13.6 are the ionization potential in 
°K and eV respectively. 

The bound-bound transitions are included for principle 
quantum numbers up to five. The absorption coefficient is 
calculated by using the relation 

= SL(v) (3.12) 

where S is the line strength and L(v) is the line shape 
factor. The line strength is given by the relation 

S = (l.lOE - 16)fn^Njj exp [ (-157780/T) (1 - l/n|) ] (3.13) 

The line shape factor is given by the relation 

L(v) = + (v - (3.14) 

where is the frequency at the line center and a is 

the line half-width, and these are given by 

= 13.6[(l/nJ) - (l/n^)l 

Y = a[1.05E 15(n^ - 

The constant a in the above equation is taken to be 0.642 
for the first line and unity for the remaining lines. 

The absorption coefficients for the free-free and bound- 
free transitions of the negative hydrogen are 

K?. = (6.02E - 39)N„N /v ^ (3.17) 

rr He 

K^f = (2.89E - 17) (3** - 43^ + 3.643^ + 0.733)Njj- (3.18) 

where 3 = 1.502/v. The threshold for the bound-free transi- 
tion of H is 0.757 eV. 


(3.15) 

(3.16) 
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The absorption coefficient for molecular hydrogen in the 
jth step is obtained in accordance with equation (3.9) and is 
expressed as 

<f = fj(T)NH^ (3.19) 

where f ^ (T) is dependent on the particular step. The molecular 
bands cover the steps from 7 to 17 ev. 

Further details on constructing the step-function model and 
utilizing it in the radiative flux equations are given in 
references 36 and 37. 

4. PHYSICAL CONDITIONS AND DATA SOURCE 

For this study, the entry body is considered to be a 45® 
hyperboloid blunt body which enters the Jovian atmosphere at a 
zero degree angle of attack. The body nose radius, R*, is 
taken to be 23 cm. The body surface is assumed to be gray having 
a surface emittance of 0.8, and the surface temperature is taken 
to be uniform at 4,000 ®K. 

4.1. Free-Stream Conditions 

Information on Jupiter's atmospheric conditions are avail- 
able in reference 45, and some of these are illustrated in 
figure 4.1. For different altitudes of entry, the free-stream 
conditions used in this study are given in reference 12 and are 
simmarized in table 4.1. The free-stream atmospheric composi- 
tion is assumed to be 85 percent hydrogen and 15 percent helium 
by mole fraction. The temperature of the atmosphere is taken 
to be constant at 145 ®K, and the free-stream enthalpy is given 
by the relation = 1.527 RT^. The number density of hydrogen 
can be calculated by the relation 

Njj^ = (7.2431172 X 10”) (p„/T^)Xjj^ (4.1) 
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where X„ is the mole fraction of Hz and p has units of 
Hz 

N/m^ . 


Table 4.1. Altitudes of entry and free-stream 
conditions for Jovian entry. 


z 

km 

V 

00 

cm/sec 

Poo 

g/cm® 

T 

00 

K 

Poo 

dync/cm^ 

p V® 
*^00 00 

£ 

116 

3.909E6 

4.65E-7 

145 

2.44E3 

2.777E13 

0.006645 

143 

4.517E6 

1.27E-7 

145 

6 . 66E2 

1.17E13 

0.01272 

190 

4.736E6 

1.33E-8 

145 

69 

1.412E12 

0.03930 

225 

4.756E6 

2.50E-9 

145 

13 

2.69E11 

0.09064 

261 

4.758E6 

4.53E-10 

145 

2.38 

4.879E10 

0.2129 


4.2. Gaseous Composition of Precursor 
and Shock-Layer Regions 

Initially the gas composition ahead of the shock is, of 
course, the free-stream composition. Absorption of high intensity 
radiation from the shock, however, produces H, Hz, and electrons 
in the precursor region as a result of photodissociation and 
photoionization. Any other species which may be produced usually 
are neglected. Thus, the precursor gas is considered to be a 
mixture of Hz» He, H, Hz^, and e“. Further information on pro- 
duction of these species in the precursor region is available 
in reference 12. 

The shock-layer gas is assumed to be in chemical equilibrium 
and is taken to be a mixture of eight chemical species, Hz, H, 

^ 4* i I 

H , H , He, H^, He , and e“. The number density of these 
species are obtained by considering five chemical reactions as 
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1. 

H2 

t 2H 


2. 

H 

t + e“ 


3. 

He 

% He"*" + e” 


4. 

He’’’ 

X He"’”’’ + e“ 


5. 

H" 

i H + e“ 

(4.2a) 

general 

, these reactions can be expressed by 


L ' 

a. A. 

JL 1 

-t- Z! ^i®i 

(4.2b) 


The number densities (particle/m^ ) are related to the equilibrium 
constant for the jth reaction by 

K. = (nN*’i(Bj^)]/[HN®i(A^)] (4.3) 

The equilibrium constants are calculated from the species partition 
function, and for the five chemical reactions they are found to 
be (ref. 36) 

)ln Ki = 52.2042 + 0.5 £n T 

+ Hn[l - exp (-6331/T) ] - 51964/T 
^,n K 2 = 35.4189 + 1.6 S,n T - 157810/T 

JJ-n Ka = 36.8052 + 1.5 2,n T - 285287/T 

Jin K 4 = 35.4189 + 1.5 J^n T - 631310/T 

Jin K = 36.8052 + 1.5 Jin T - 8750/T (4.4) 

The conservation equations for the hydrogen and helium 
nuclei and charge are 

•’h + * ”h- = “h 
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(4.6) 




He 


^He+ ■*■ “ ^H“ ^e“ 


(4.7) 


The number densities of the hydrogen and helium nuclei are 
calculated by 




(4.8) 


^He ^He^^o^'^^o^ 


(4.9) 


where 


M = 2.016x„ + 4.003x„^ 

o H 2 He 

In the above equations, represents the Avogadro's constant, 

p is the mixture density in g/cm^ , x „2 is the mole fraction 

of molecular hydrogen, and x„ is the mole fraction of helium. 

The solution procedure for obtaining the eight unknown 
number densities is discussed in reference 36. The closed-form 
solutions are obtained by solving equation (4.3) for each reaction 
independently. This is accomplished by setting the appropriate 
values in equations (4.5) to (4.7) to zero if the species are 
not present in the reaction. The closed-form solutions for the 
number densities (in particles/cm^ ) of each species are given 
by 


H 2 ; 

^H2 

= (N°/2) + (Ki/8)[(1 + 8N°/Ki)-^'^^ - 1] 

(4.10) 

H-*-: 

^H+ 

= (K/2)[(l + 4N°/K2)^'^^ - 11 

(4.11) 

H: 


= N° - 

(4.12) 


28 



1] 




Njj+ = (Di/ 2 ) [(1 + 4K3NH3/D!) 


1/2 _ 


Di = K3 + Njj+ 


(4.13) 


+ + 
He 


N 


He 


(Dz/2)[(1 + 4 K^N°^/d|)^/^ - Hr 

D2 = K4 + Njj + 


(4.14) 


He; 


N = N - N + N ++ 
^He ^He He He 


(4.15) 


N _ 

e 




(4.16) 


H : N„- = N„N -/Ks (4.17) 

H He 

4 . 3 Thermodynamic and Transport Properties 

Thermodynamic properties for specific heat, enthalpy, and 
free energy, and transport properties for viscosity and thermal 
conductivity are required for each species considered in different 
flow regimes. For the precursor zone as well as shock layer, 
the general expression for total enthalpy, specific enthalpy, 
and specific heat at constant pressure are given respectively by 


= h + (u^ + v^)/2 

(4.18) 

■ ^ X . h . — ^ ' (C /M ) h 

(4.19) 

= E x.C . 
^ 1 pi 

(4.20) 


However, specific relations for h and 
are quite different. 


for the two regions 


For the precursor region, the relation for the specific 
enthalpy is obtained by following the procedure described by 
Smith (ref. 1) as 
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h = 1.4575RT + (0.75RT + D)C^ + (1.25RT + 1/2)C^^ 


(4.21) 


where D and I represent the dissociation and ionization 
energy respectively, and their values are available in 
reference 40. 

The derivation of equation (4.21) essentially follows from 
the consideration of equation (4,19). If it is assumed that the 
internal energy of each particle can be described only by 
translational and rotational modes, then the relation for specific 
enthalpy of each species can be expressed as 

^He = I RT + p/p = I RT 

= •§ + I RT + p/p + I = |- RT 

hjj+ = I RT + I RT + p/p + I = ~ RT + I 

3 5 

hfl = 2 p/p + D = ^ RT + D 



Also, from the conservation of charged particles one can 
write 


Now , for 85 percent H 2 and 15 percent He on volume basis (or 
76 percent H 2 and 24 percent He on mass basis), equation (4,19) 
is written as 


30 



X:(cyM^)h^= (0.26/4) (5RT/2) + [(0.74 - C^+ - Cjj)/2] (7RT/2) 


+ [(5RT/2 + D)]Cjj + (7RT/2 + 1) {C^+/2) 

+ (5RT/2) (Cjj+/2) 

A simplification of the above equation results in equation (4.21). 

In the shock region, equations (4.19) and (4.20) are used 
to calculate h and C^. With representing the mole 

fraction of the ith species, the expressions for h^ and 
are found from references 46 and 47 as 

h^ = RT[ai + (a2/2)T + (a3/3)T^ + (a4/4)T^ 

+ (a5/5)T** + ae/T] (4.22) 

C . = R(ai + a 2 T + agT^ + a^T^ + asT*" (4.23) 

P ^ 

where R is the universal gas constant (= 1.98726 cal/mole - K) 
and T is the local fluid temperature in K. For different 
species, values of the constants ai, a 2 / ... are given in 

reference 47, and for species under present investigation they are 
listed in table 4.2. It should be pointed out here that in this 
study, instead of employing equation (2.17b), equations (4.18), 
(4.19), and (4.22) are used to calculate the enthalpy variation 
in the shock layer. This is because slightly better results are 
obtained by using the above set of equations. 

For the shock-layer gas, the mixture viscosity and thermal 
conductivity are obtained by using the semi-empirical formulas 
of Wilke (ref. 48) as 
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u> 

M 


Table 4.2 


Coefficient for evaluation of the specific heat at constant 
pressure and enthalpy for various hydrogen/helium species. 


Species 

ai 

^2 

Coefficients 
a 3 a 1 , 

^5 

ae 


Temp. 

Range 

K 


2.5 

0 

0 

0 

0 

2.547162E+4 

> 

300 

H 

2.5 

0 

0 

0 

0 

2.547162E+4 

> 

1,000 


2.475164 

7.366387E-5 

-2.537593E-8 

2.386674E-12 

4 .551431E-L7 

2.523626E+4 

> 

6,000 


3.057445 

2.676520E-3 

-5.809916E-6 

5.521039E-9 

-1.812273E-12 

-9.889047E+2 

> 

300 

Hz 

3.10019 

5.111946E-4 

5.264421E-8 

-3.490997E-11 

3.694534E-15 

-8.773804E+2 

> 

1,000 


3.363 

4.656000E-4 

-5.127000E-8 

2.802000E-12 

4.905000E-17 

-1.018000E+3 

> 

6,000 


2.5 

0 

0 

0 

0 

1.840334E+5 

> 

300 

h'^ 

2.5 

0 

0 

0 

0 

1.840334E+5 

> 

1,000 


2.5 

0 

0 

0 

0 

1.840334E+5 

> 

6,000 


2.5 

0 

0 

0 

0 

-7.453749E+2 

> 

300 

He 

2.5 

0 

0 

0 

0 

-7.453749E+2 

> 

1,000 


2.5 

0 

0 

0 

0 

-7.453749E+2 

> 

6,000 


2.5 

0 

0 

0 

0 

2.853426E+5 

> 

300 

He 

2.5 

0 

0 

0 

0 

2.853426E+5 

> 

1,000 


2.5 

0 

0 

0 

0 

2,853426E+5 

> 

6,000 


2.5 

0 

0 

0 

0 

-7.453749E+2 

> 

300 

e 

2. 5 

0 

0 

0 

0 

-7.453749E+2 

> 

1,000 


2.508 

-6.332000E-6 

1.364000E-9 

-1. 094000E-13 

2.934000E-18 

-7.450000E+2 

> 

6,000 



(4.24) 


N N 

i=l ^ j=l 


X . (}) . . ) ] 

D ID 


N N 

K=Etx.K./(J^ x.())..)] (4.25) 

i=l ^ ^ j=l 3 13 

where 

. = [1 + (M /M )^/'^]V{VS[i + (M./M 

-1»J -*-J J1 XJ 

and is the molecular weight of species i. For hydrogen/ 

helium species, specific relations for viscosity and thermal 
conductivity are given in references 49 and 50. The viscosity 
of Ha and He, as a function of temperature, can be obtained from 
reference 49 as 

y„ = (0.66 X 10“M (T) ^^^/(T + 70.5) (4.26) 

tl2 

y„ = (1.55 X 10“® ) (T) ^/^/(T + 97.8) (4.27) 

He 

The thermal conductivity of Ha and H are obtained from reference 
50 as 


K„ = 3.212 X 10-^ + (5.344 x 10"®)T (4.28) 

Ha 

K„ = 2.496 X 10-^ + (5.129 x 10"®)T (4.29) 

II 

The viscosity of H and thermal conductivity of He are obtained 
from the relation between viscosity and thermal conductivity 
of monatomic gases as given in reference 48 by 

K = (15/4) (R/M)y (4.30) 

Very little information is available on transport properties of 
other species such as Ha^, H^, e“ , etc. Fortunately, transport 


33 



properties are important only in the boundary-layer region where 
the temperature is not high enough to produce these species. 

It should be noted that all relations presented in this 
section are expressed in dimensional form. 

5 . SOLUTION PROCEDURES 

An iterative procedure has been used to couple the precursor 
and shock-layer solutions. In this method, the shock-layer 
solutions are obtained first with no consideration of precursor 
effect. From this solution, the radiative flux at the shock front 
(which influences the precursor region flow) is determined. By 
employing this value of the radiative flux, different precursor 
region variables are calculated through use of equations (2.6) 
through (2.10). Values of these flow variables are obtained 
just ahead of the bow shock, and then the Rankine-Hugoniot 
relations are used to determine the conditions behind the shock. 
These conditions are used to obtain new shock-layer solutions 
from which a new value of the radiative flux at the shock is 
calculated. The procedure is continued until the radiative 
flux at the shock becomes invariant. 

The solution procedures for the precursor and shock-layer 
regions are described in some detail in the following subsections. 

5.1. Precursor Region Solutions 

A direct integration of equations (2.6) through (2.10) 
results in the following governing equations for the precursor 
region 


pv = 


p V (u - u ) = 0 


Poo ’^oo^v - + (p - P^) = 0 


(5.1) 

(5.2) 

(5.3) 
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0 


(5.4) 


CO 00 00 R 

P v_(9C^/3n) - K = 0 

*^0O 00 ' Qt' ' Qt 


(5.5) 


In equation (5.4), H represents the total enthalpy and is 
given by a combination of equations (4.18) and (4.21). The 
expression for the radiative flux, q^, is given by equation 
(3.4). For the present application, equation (5.5) will be 
written for atomic hydrogen and hydrogen ions . By following 
the procedure outlined in reference 12 and 19 the expressions 
for species concentration are found to be 


N 

= -264 E Y 
i=l 


D_E3(Ti)(kTs) 


-1 


I(vJ) 


(5.6) 


N 

C„+ = - 264EY-T E3(tJ (kT^) - (v?) (5.7) 

H 2 1 s 1 


where 


64 = (15/7t'») [q(0)mi/(p^V^)] 


I (v?) 



V 


li 


{v^/ [exp (v) 


1] }dv 


and mx represents the weight of the Ha molecule in grams per 
molecule . 

Note that there are nine algebraic equations to evaluate the 

nine unknowns p, v, u, p, h, H, C„+* The solutions 

H hia 

of this set of algebraic equations are obtained by using the 
Gauss-Seidel method (ref. 51). The properties at the infinity 
are used as the first initial guess in the Gauss-Seidel method. 
The iteration is continued until all the quantities in this 
region become invariant. The flow chart of the computational 
procedure is illustrated in figures 5.1 and 5.2. 
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5.2. Shock-Layer Solutions 


A numerical procedure for solving the viscous shock-layer 
equations for stagnation and downstream regions is given by 
Davis (ref. 22). Moss applied this method of solution to 
reacting multicomponent mixtures in references 9 and 10. A 
modified form of this procedure is used in this study to obtain 
solutions of the viscous shock-layer equations. In this method, 
a transformation is applied to the viscous shock-layer equations 
in order to simplify the numerical computations. In this trans- 
formation most of the variables are normalized with their local 
shock values. The transformed variables are (ref. 9): 


n = y/Yg 

p = P/P 3 

P = P/Pg 


S = X 

CO 

CL 

a 

II 

1 Q. 

K = K/K 

S 


3 = U/Ug 

T = T/Tg 

C = C /C ^ 
P P PS 


3 = v/v^ 

H = H/H 

s 


(5.8) 


The transformations relating ' the differential quantities are 


— ( ) 
9x ^ ^ 


(dy^/aa.^ ( > 


(5.9a) 


where 


( ) 

ay ^ ^ an ^ 


— ()=— — 


{ ) 


(5.9b) 


ay' 




The transformed equations can be expressed in a general form as 


a^w/an^ t a i aw/ an t s-2W + as + a 4 aw/ a ^ — 0 


(5.10) 


The quantity W represents u in the x-momentum equation, T in 
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the temperature energy equation, H in the enthalpy energy equa- 
tion, and in the species continuity equations. In most cases, 

the coefficients ai to a 4 to be used in this study are exactly 
the same as given in reference 9. However, there is one exception. 
Since radiation is included in the present study, the coefficients 
of the energy equation are different from those used in reference 
9. For example, in the enthalpy energy equation, coefficients 
aj , a 2 / and a 4 are the same as given in reference 9, but 
as is different, and this is given by 


a 3 


P P y 
r , s r-^ s 


P yH 
s s 




( _!£ 

+ YgHK 
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cos 0 

+ y n cos 9 
s 
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y P p V V 
•^s r s s 




lE 

3n 


Us 


where 


ip = 
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,sL p 


cos 0 


YghK 


r + 




cos 0 


)] 




i=l 


dC. 
1 

3n 


+ 


u^yu 

Pr 


(Pr,^P 
' s r 


- 1 ) 


ini 

9nJ 


(5.11) 


1 + YgHic 


Other transformed equations are the same as given in reference 9. 

The surface boundary conditions in terms of transformed 
variables are 


y = 0, v=0, T = T 

w 

The transformed shock conditions are found to be 


(5.12) 


U = v= T = H= p = p = l 


(5.13) 


at n 


1 . 
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The second order partial differential equations as expressed 
by equation (5.10), along with the surface boundary and shock 
conditions, are solved by employing an implicit finite-difference 
method. In order to obtain numerical solutions for the down- 
stream region, it is necessary to have an accurate stagnation 
streamline solution. Since the shock shape is affected by the 
downstream flow, a truncated series of shock standoff distance 
is used to develop the stagnation streamline equations. As such, 
the shock standoff distance is expressed by 

= yis + 

Since ^ is small and the curvature k is approximately one 
in the stagnation region, it is logical to say that (see fig. 2.1) 

3 ^ (5.15) 

Since 0 = ('ir/2) - 3, there is obtained 

ot = 0 + tan“^ [ On^/a?) /(I + Ky^) ] 

= (it/ 2) + 4{[2y2g/(l + Vig)] - l} (5.16) 

By using equations (5.14) to (5.16), the shock relations 
[eqs. (2.21) - (2.26)] can be expressed in terms of 


expanded variables as 



(5.17) 

u*_ = -5[2y^^/(l + y^g) - 1] 

(5.18) 

“s- = 5{l - (2y,3/(l + y: 3 )](l + I/P 3 -)) 

(5.19) 


Pg- = Pg+ + (1 - 1/Pg-) + (1 - 1/Pg-) 

• u - 2Y^g/(l + y,s)l"} 

V = ’’s+ + (1 - l/Pa-)/2 
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(5.21) 



In equations (5.17) through (5.21), only p and u involve 
in the first terms of their expansions. Thus, a series 
expansion for the flow variables is assumed about the axis of 
symmetry with respect to nondimens ional distance ^ near the 
stagnation streamline as 


P(^/H) 

= Pi (ti) p2 (n) + 

(5.22a) 

u(S,ri) 

= uj (n) C + 

(5.22b) 

<! 

= v(n) + 

(5.22c) 

piKfT)) 

= Pi (n) + 

(5.22d) 

T(^,n) 

= Ti (n) + 

(5.22e) 

y (Cf n) 

= y 1 (n) + — 

(5.22f) 

K(^,n) 

= Ki (n) + — 

(5.22g) 

Cp ( ^ , n ) 

= c (n) + 

pi 

(5. 22h) 

c. ( 5 ,n) 

= c. (n) + — 

(5.22i) 


Since y„ is a function of downstream flow, it cannot be deter- 
■‘2 s 

mined by the stagnation solutions. Thus, a value of y^^ = 0 
is assumed initially. This assumption is removed by iterating on 
the solution by using the previous shock standoff distances to 
define Y 2 s* 

The new form of x~momentum and energy equations in the ^,n 
plane can be written as 


— — — 


ai 


aw 

an 


+ a2W + as 


0 


(5.23) 


For x-momentum, W = u and coefficients ai, a. 2 , and as are 
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exactly the same as given in reference 9. For the enthalpy 
equation, W = H and again ai and a 2 are the same as defined 
in reference 9, but as is given by 


+ 2 Y /(1 + 




(Vi 


P p V Vi/e 
s r^s s 


77 2 , . — ■ 


y^pH^) (9pi/3n) 


(5.24) 


Other stagnation streamline equations are the same as given in 
reference 9. It should be noted that at the body surface 
Pi = 1 and p2 =0. 

As mentioned earlier, the governing second-order partial 
differential equations are solved by employing an implicit finite- 
difference method. For this purpose, the shock layer is considered 
as a network of nodal points with a variable grid space in the 
rj-direction . The scheme is shown in figure 5.3, where m is a 
station measured along the body surface and n denotes the 
station normal to the body surface. The derivatives are converted 
to finite-difference form by using Taylor's series expansions. 

Thus, unequal space central difference equations in the n-direction 
at point m, n can be written as 
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(5.25b) 
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I 


w - w 

3W » _ m/ii m-ljTi 


(5.25c) 


A typical difference equation is obtained by substituting the 
above equations in equation (5.10) or (5.23) as 


W 


m,n 


- (D /B ) - (A /B )W - (C /B )W „ T 

n n n n m,n+l n n iri/n-1 


(5.26) 


where 


A^ = (2 + aiAn^_^)/[Anj^(Anj^ + 


= - [2 - ai(in„ - - a2 - 


c„ = (2 - aiAn^)/[An„_i(An„ + an„_i)] 


= a, - 


Now, if it is assumed that 


r^ = E W . + F 

HI / T1 HI / n * ^ J1 


(5.27) 


or 


W , = E tW + F , 
m,n-l n-1 m,n n-1 


(5.28) 


then by substituting equation (5.28) into equation (5.26), there 
is obtained 


W. 


m 


,n = 


+ (-D - C F ,)/(B + C E ,) 

' n n n-l^'^ n n n-1' 


(5.29) 
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By comparing equations (5,27) and (5.29), one finds 


E =-A/(B +CE t) 
n n n n n-1 


(5.30) 


F 

n 


(- - 


C F 1 )/(B + C E , ) 

n n-1 ' n n n-1 


(5.31) 


Now, since Ei and Fi are known from the boundary conditions, 

E^ and F^ can be calculated from equations (5.30) and (5.31). 

The quantities W at point m, n can now be calculated from 

m,n ^ 

equation (5.27). 

The overall solution procedure starts with evaluation of the 
flow properties immediately behind the shock by using the 
Rankine-Hogoniot relations. With known shock and body surface 
conditions, each of the second-order partial differential equations 
are integrated numerically by using the tridiagonal formalism of 
equation (5.10) and following the procedure described by equations 
(5.26) to (5.31). As mentioned before, the solutions are obtained 
first for the stagnation streamline. With this solution providing 
the initial conditions, the solution is marched downstream to the 
desired body location. The first solution pass provides only 
an approximate flow-field solution. This is because in the first 
solution pass the thin shock-layer form of momentum equation is 
used, the stagnation streamline solution is assumed to be inde- 
pendent of downstream influence, the term dy^/d? is equated to 
zero at each body station, and the shock angle a is assumed to 
be the same as the body angle 6. All these assumptions are 
removed by making additional solution passes. 

In the first solution pass, the viscous shock-layer equations 
are solved at any location m after obtaining the shock conditions 
(to establish the outer boundary conditions) from the precursor 
region solutions. The converged solutions at station m-1 are 
used as the initial guess for the solutions at station m. The 
solution is then iterated locally until convergence is achieved. 
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For the stagnation streamline, guess values for dependent 
variables are used to start the solution. In the first local 
iteration, both and OW/BC) are assumed to be zero. 

The energy equation then is integrated numerically to obtain 
a new temperature. By using this temperature, new values of 
thermodynamic and transport properties are calculated. Next, 
the x-momentum equation is integrated to find the u component 
of velocity. The continuity equation is used to obtain both the 
shock standoff distance and the v component of velocity. The 
pressure p is determined by integrating the normal momentum 
equation. The equation of state is used to determine the density. 

For example, the integration of the stagnation streamline 
continuity equation from 0 to n results in 

[(1 + y,sn)^P,3V^^p,]v, = (-2y,3P,sU,g)A (5.32) 


where 



+ y igTi) Piuidn 


(5.32) 


This equation gives the v-velocity component along the stagnation 
streamline. However, integration of the continuity equation from 
n = 0 to n = 1 results in 


(1 + 


^is^ 


P V 

1 S IS 


= - 2p u, y, (B 
IS ISIS 


+ C) 


(5.33) 


where 

1 1 

B = J' piudn, c = y f piuindn 
o o 

The shock standoff distance can be obtained from the solution of 
equation (5.33) as 
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(5,34) 


^ -(2v,^ ZBU'^l -I- C(2v,^ - 4(v,^ + 

Similarly, other quantities at the stagnation streamline are 
obtained. 

With known stagnation streamline solution and body surface 
and shock conditions, the above procedure is used to find solutions 
for any body location m. The downstream shock standoff distance 
and the v-velocity component are obtained by integrating the con- 
tinuity equation in the n-direction from 0 to 1, and 0 to n 
respectively. Integration of the continuity equation from 
n = 0 to n = 1 results in 

1 1 

■^[y| cos e PgU^ f pundn + y^rp^Ug f pudn 


= (r + y^ cos 6 ) [ygPg^s “ (1 + (5.35) 

By defining, for station m 

Cl = cos 6 Pg^g ^ puridri , Cz 

and denoting the same relations by C 3 and C 4 ^or station 
m- 1 , equation ( 5 . 35 ) can be expressed in terms of a difference 
equation as 


= rp u^ 
•^s s 


J 
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pudn 


[(C.y^ + - (C3y^ + C,yg)„.iJ/i5 




- rp V Ky - cos 9 P^v^y^^ - cos 6 p^v Ky 
K ^ a -I ejTTi ^ S 3 -^ sm <? e; 


s s sm 


s s sm 


(5.36) 
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This can be expressed in a quadratic form as 


sm sm 


= 0 


(5.37) 


where 


AA = Cl + cos 9KpgVgA^ 


EB - Cz + rpgVgKAC - cos 0 PgU^y^A^ + cos 0 PgV^A^ 


CC = - 




The shock standoff distance at station m is obtained from 
equation (5,37) as 


y^™ = {-(BB) + [(BB)^ - 4(AA) (CC)l^/^}/2(AA) 


sm 


(5.38) 


The v-velocity component can be obtained in a similar manner. 
Integration of the continuity equation from 0 to y gives 


-^r f y (r + y n cos 0)p u pudr;! 
9 ? sm sm ^ s s J 


+ (r + y^^n cos 0) [(1 + 

■ ysm’1Ps“sP"' “ ° 

As before, this can be expressed in terms of a difference 
equation as 

{[(KK) - (KK)„_J/A5) + (II)„V + (JJ)^ 0 (5.40) 

where 


45 



(II) = (r + y n cos 0) (1 + y„^nK)p^v p 
m sm sm s s 

(JJ) = “ (r + y cos 6)y' np u pu 
' -^sm ' ^-^sm '^s s 

n 

(KK) = / y„_(r + y ri cos 0)p u„pudr| 
m ~'-^sm -^sm' ^ss*^ ' 

o 

ThuSr the v-velocity component at each point on the station 
m can be obtained from equation (5.40). Other quantities at 
station m are obtained by a similar manner. As mentioned 
before, the first pass is only an approximate solution because 
of several inherent assumptions. These assumptions are removed 
by iteration in the next pass. For the subsequent solution 
passes, the shock angle and y^^ are given by 

a = 0 + tan”' ] (5.41) 

^32 " <^33 " ys,)/4(AO" (5.42) 

The flow diagrams for computation procedure are shown in figures 
5.1, figures 5.4 to 5.8. 


6. RESULTS AND DISCUSSION 

The governing equations of both the precursor and shock-layer 
regions were solved for physically realistic Jovian entry condi- 
tions. Results of the complete parametric study are presented 
in this section. First, the results are presented for quantities 
just behind the shock wave, and then a few results of flow 
variables within the shock layer are presented. Next, results 
are presented for the entire shock-precursor region. Finally, 
a few results are presented to demonstrate the influence of 
precursor heating (i.e., preheating of the gas in front of the 
shock) on different heat fluxes. 

The radiative flux from the shock layer towards the precursor 
region is found to be highest at the stagnation line shock location. 
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Results of the radiative flux from the shock front are shown in 
figure 6.1 for different altitudes of entry. As would be expected, 
precursor heating results in a higher radiative flux at the shock 
front. It is seen that the radiative flux reaches a maximum value 
for an altitude of about 116 km, and the largest precursor effect 
(PE) of about 8 percent is found to be for this altitude. This 
is a direct consequence of the free stream and entry conditions 
at this altitude. For other entry conditions (altitudes) , 
precursor effects are seen to be relatively lower. 

Figure 6.2 shows the shock standoff variation with distance 
along the body surface for different entry altitudes. The shock 
standoff distance, in general, is seen to decrease with increasing 
altitudes. This is because higher entry velocities are associated 
with higher altitudes. The precursor heating results in a slight 
increase in the shock standoff distance (a maximum of about 2 
percent for z = 116 km) because the density of the shock layer is 
slightly reduced. 

The conditions just behind the shock are illustrated in 
figures 6.3 to 6.8 as a function of distance along the body for 
different entry altitudes. For z = 116 km, figure 6.3 shows that 
precursor heating increases the enthalpy by a maximum of about 
2 percent at the stagnation streamline. The change in shock 
temperature is shown in figure 6.4 for different altitudes. As 
would be expected, precursor heating results in a relatively 
higher temperature. The effect of precursor heating on the 
pressure just behind the shock is shown in figure 6.5, and it 
is seen to be very small. Since the pressure essentially remains 
unchanged, precursor heating results in a decrease in the density, 
(see fig. 6.6) mainly because of an increase in the temperature. 
Figure 6.7 shows that precursor heating has no significant 
influence on the u-velocity component, but the v-velocity is 
slightly increased (see fig. 6.8) as a result of decrease in 
the shock density. 

Variations in temperature, pressure, density, velocity, chem- 
ical species, and transport properties across the shock layer 
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are shown in figures 6.9 to 6.13 for entry conditions corresponding 
to an altitude of z = 116 km. Results presented in these figures 
are normalized by their shock values, and they show that pre- 
cursor effects are felt throughout the shock layer. Results 
presented in figures 6.9 to 6.11 for two body locations (? = 0 
and 1) indicate the relative change in temperature, pressure, 
density, and velocity as compared with their shock values. 

For ^ = 0, figure 6.12 shows that precursor heating slightly 
decreases the concentration (mole fraction) of atomic hydrogen 
and increases the concentration of ions and electrons throughout 
the shock layer. For the stagnation streamline, results pre- 
sented in figure 6.13 indicate that the transport properties 
(y and K) are changed only slightly in the vicinity of the 
bow shock. 

Variations of temperature, pressure, density, and velocity 
along the stagnation streamline in the entire shock-layer/precursor 
zone are illustrated in figures 6.14 to 6.17 for different 
altitudes. Since higher entry velocities are associated with 
higher altitudes, precursor effects, in general, are found to 
be larger for higher altitudes. The results for the precursor 
region show a dramatic increase in the pressure and temperature 
but only a slight change in the density and velocity. The 
changes are largest near the shock front because a major portion 
of radiation from the shock layer gets absorbed in the immediate 
vicinity of the shock front. A complete discussion on changes 
in flow variables in the precursor region is given by Tiwari 
and Szema in reference 12. Figures 6.14 and 6.15 show that, 
in spite of a large increase in the temperature and pressure 
in the precursor region, precursor heating does not change 
the temperature and pressure distribution in the shock layer 
dramatically. The change in temperature, however, is signifi- 
cant, the maximum change occurring, as would be expected, just 
behind the shock. There is a slight change in the pressure 
near the body, but virtually no change closer to the shock. 

Figure 6.16 shows that the change in density in the shock 
layer is higher for higher altitudes and towards the shock. 
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As discussed before, precursor heating results in a slight decrease 
in the shock- layer density. Virtually no change in the u-component 
of the shock-layer velocity was found, but, as shown in figure 
6.17, the v-component is slightly increased. 

The effects of precursor heating on different heat fluxes 
in the shock layer are illustrated in figures 6.18 to 6.24. 

These results clearly demonstrate that precursor heating has a 
significant influence on increasing the heat transfer to the 
entry body. This increase essentially is a direct consequence 
of higher shock- layer temperatures resulting from the upstream 
absorption of radiation. Figures 6.18 to 6.20 show the variations 
in radiative and convective heat fluxes with distance along the 
body surface for three different altitudes. For z = 116 km, 
results illustrated in figure 6.19 reveal that the precursor 
heating results in a 7.5 percent increase in the radiative flux 
and about a 3 percent increase in the convective flux to the 
body at the stagnation point. The increase in heat transfer at 
other body locations and for other entry conditions (altitudes) 
is relatively lower. A similar conclusion can be drawn from the 
results presented in figures 6.21 to 6.23 for the radiative flux 
towards the shock and the body for two body locations (^ = 0 and 
1) and for three different entry conditions. Results of radiative 
and convective heat flux at the body (for ^ = 0) are illustrated 
in figure 6.24 for different altitudes of entry. The radiative 
flux results are seen to follow the trend exhibited in figure 
6.1 for radiation at the shock front. The convective heat flux, 
however, is seen to increase slowly with the altitude up to 
z = 131 km, and thereafter to decrease with increasing altitudes. 
The precursor effect is found to increase the radiative heating 
by a maximum of about 7.5 percent at z = 116 km and the convec- 
tive heating by 4.5 percent at z = 131 km. 

It should be emphasized here that the results presented in 
this study were obtained for carefully selected entry conditions. 
For other conditions, precursor heating may have an entirely 
different influence on flow variables in the shock layer. For 
example, if entry velocities higher than 38 km/s are considered 
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for altitudes lower than 116 km, then one would expect precursor 
heating to have a considerably greater influence on the shock 
layer flow phenomena. This, in turn, would result in higher 
radiative and convective heating of the entry body. 

7. CONCLUSIONS 

The main objective of this study was to investigate the 
influence of precursor heating on the entire shock-layer flow 
phenomena. For Jovian entry conditions, results were obtained 
for the precursor zone as well as the shock layer. For the pre- 
cursor region, results indicate that the temperature and pressure 
are increased significantly because of absorption of radiation 
from the shock layer, but only a slight change is noticed in 
other flow variables. In the shock layer, results of flow 
variables were obtained along the body and the bow shock and 
across the shock layer for different altitudes of entry. The 
influence of precursor heating was found to be larger at higher 
altitudes. Specific results indicate that, in most cases, pre- 
cursor heating has a maximum influence on flow variables (except 
the pressure) at the stagnation line shock location. It was 
found that while pressure essentially remains unchanged in the 
shock layer, the precursor heating results in an increase in 
the enthalpy, temperature, and v-component of velocity, and a 
decrease in the shock layer density. For the entry conditions 
considered in this study, results clearly demonstrate that pre- 
cursor heating has a significant influence on increasing the 
heat transfer to the entry body. It was found that the radiative 
heating is increased by 7.5 percent at z = 116 km and the convec- 
tive heating by 4.5 percent at z = 131 km. 

For further study, it is suggested that the precursor region 
flow phenomena be investigated without making the thin layer 
approximation. Also, it would be advisable to use different 
free-stream temperatures for different entry altitudes. In the 
shock layer, the case of chemically nonequilibrium flow should 
be included, and then the influence of precursor heating on 
different flow variables should be investigated. 
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ABSORPTION CROSS-SECTION 








PRECURSOR SHOCK-LAYER SOLUTION PROCEDURE 



Figure 5.^1 


Flow chart for combined precursor/ 
shock-layer solution procedure. 

















CTl 

K) 



Figure 5.3. Finite difference representation of flow field. 




SUBROUTINE SHOCK 


Figure 5.4. Flow chart for subroutine SHOCK for 
shock-layer solution. 
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SUBROUTINE SHOKLY 



Figure 5.5. Flow chart for subroutine SHOKLY fob shock- 
layer solution. 
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SUBROUTINE ENERGY 



Figure 5.6. Flow chart for subroutine ENERGY for 
shock-layer solution. 
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SUBROUTINE MOMENTM 



Figure 5.7. Flow chart for subroutine MOMENTM for 
shock-layer solution. 


66 
















SUBROUTINE 



subroutine RADIATION for 
ution. 
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INTEGRALS FOR SUBROUTINE RADIATION 



Figure 5.8b. Definition of integrals used in subroutine 
RADIATION. 
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Figure 6.1. Radiation flux towards the precursor region 
at the stagnation line shock location. 





Figure 6.3. Enthalpy variation just behind the shock with 
distance along the body surface. 









— PRECURSOR EFFECT 

— NO PRECURSOR EFFECT 



Figure 6.7. Variation of u- velocity component just behind 
the shock with distance along the body surface 







Figure 6.8. Variation of v-velocity component just behind 
the shock with distance along the body surface 


T*/T 



n = y/Yg = y*/yg 

Figure 6.9. Variation of temperature in the shock layer for two body locations - 0 and 1) 
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Figure 6.11. Variation of velocity components in 
locations (^ = 0 and 1} . 





PRECURSOR EFFECT 


NO PRECURSOR EFFECT 





Figure 6.14. Temperature variation in the shock/precursor 
region along the stagnation streamline. 



Figure 6.15. Pressure variation in the shock/precursor 
region along the stagnation streamline. 









p = P*/P5, 
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Figure 6.16. Density variation in the shock/precursor 
region along the stagnation streamline. 



Figure 6.17. Variation of v-velocity component in the 

shock/precursor region along the stagnation 
streamline . 













Figure 6.22. Variation of radiative heat flux in the shock 
layer for two body locations (C = 0 and 0.5), 
z = 116 km. 
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APPENDIX A 


EXPLANATION OF SYMBOLS USED IN COMPUTER 
PROGRAM OF THE SHOCK-LAYER REGION 



APPENDIX A 


AREA 

CH 

CHE 

CHP 

CHZ 

CH2P 

CE 

CPU 

CPL 

CPS 

D 

DD 

DENI 

DI 

DS 

DX 

DY 

EPS 

ETHC 

E2 

E3 

FI 


EXPLANATION OF SYMBOLS USED IN COMPUTER 
PROGRAM OF THE SHOCK-LAYER REGION 

shock angle defined in figure 2.1 
mass fraction of H 
mass fraction of He 

“ih 

mass fraction of H 
mass fraction of Hz 
mass fraction of H^ 
mass fraction of e” 
free-stream specific heat 

specific heat just ahead of shock, nondimensonal 
specific heat just behind the shock, nondimensional 
density, nondimensional 

coordinate measured normal to the body, cm 

free-stream density, cm/sec 

density just ahead of the shock, g/cm^ 

density just behind the shock, nondimensional 

step distance between two nodal points along the body 
surface 

step distance between two nodal points normal to the 
body surface 

Reynolds number parameter, nondimensional 
enthalpy, erg/g 

second exponential integral function 
third exponential integral function 
defined function of temperature 
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H2 

H 

HPLUS 

HE 

HE2PLUS 

M 

MAC I 
MV 
MVS 
MUREF 

NOPER 

NORAD 

NS 

P 

PERI 

PINF 

PS 

QO 

QFXM 

QFXP 

QTOTAL 

RAT 

T 

TB 

TH 


number density of Hz, number/m^ 
number density of H, number/m ^ 
number density of H^, number/m^ 
number density of He, number/m® 
number density of He^, number/m^ 
molecular weight 
Mach number 
viscosity 

viscosity at shock 
refrance vicosity 

! 1 without precursor effect 
2 with precursor effect 

! 1 without radiation 
2 with radiation 
shock stand-off distance 
pressure, nondimensional 
free-stream pressure 
pressure just ahead of the shock 
pressure just behind the shock 

radiative heat flux which come out from the shock 
toward the precursor region 

radiative heat flux toward body 

radiative heat flux toward shock 

net radiative heat flux between body and shock 

the ratios of dy(N + l)/dy(N) 

temperature, nondimensional 

body temperature, °K 

enthalpy, nondimensional 



THS 

TINF 

TK 

TS 

V 

VS 

V 

VEL<I)I 

VS 

VT 

XH 

XH2 

XHE 

XH2I 

XHEP 

XHE2P 

V 


enthalpy at the shock, nondimensional 
free- stream temperature, 

thermal conductivity of mixture, nondimensional 

temperature at the shock, nondimensional 

velocity component tangent to the body surface, 
nondimensional 

shock velocity ocmponent tangent to the body surface, 
nondimens ional 

velocity component normal to the body surface, 
nondimensional 

free-stream velocity 

shock velocity component normal to the body surface, 
nondimensional 


velocity at precursor region 


mole fraction of ( 


H 

between 

shock 

and 

body 

H2 

between 

shock 

and 

body 

He 

between 

shock 

and 

body 

Hz 

at free 

stream 


He"^ 

between 

shock 

and 

body 

Het 

between 

shock 

and 

body 


normal coordinate, nondimensional 


YH 

YK 

ZTA 


Planck constant, erg/s 
Boltzmann constant, erg/°K 
body angle 


96 



APPENDIX B 


COMPLETE COMPUTER PROGRAM 


PROGRAM COMPLET M NPUT » OUTPUT « TAP£5= I NPUT * TAPE6=0UTPUT ) 

COMMON /CONST 1/ AL i AK . CT , K ♦ K 1 I i MBP * MA 1 , M A2 

COMMON /PERC/ V I ( 20 ) , P I NF < 20 > * ETHP(20), D I ( 20 ) * CPL ( 20 ) . Q0l(20)* 

1 PINFN(20) 

COMMON /INFIX T I NF » PREI« VELOI» DENI* J* ALT* I 111* CPII*XH2I 
1 * NORAD *NOPRE 

COMMON /RAD/ COFl (58*51)* QFXP ( 5 1 ) * QFXM ( 5 1 ) * QFX ( 51 ) • DQY ( 5 1 ) 

COMMON /SHOCK/ OS (20)* TSC20)* VS(20)* US (20)* PS ( 20 ) * M'JS < 20 ) • TKS ( 
120 ) * ARFA (20 ) ,CPS (20 ) *THS(20) 

COMMON /SDFDIS/ NS ( 20 ) * NS 1 * NS2 * PS 1 * AS ( 20 ) *NSS2 
DIMENSION 00 (20 )* SPEC (6 ) 
real NS 

DATA (SPEC( I ) * 1 = 1 *6) / - 1 ♦ - 1 * -1 * 0 . - 1 * - 1 / 

CALL SYSTEMC (115*S*^EC) 

C 

c XH2I — MOLE FRACTION OF H 2 

C *** NORAD-1 — NO radiation **** 

c NORAr-2 — with radiation **** 

C **'< NOPRcsl — NO precursor EFFECT *#* 

C *** N0PRE=2 — WITH precursor EFFECT *** 

C 

N0PRE=2 
NORAD=2 
XH2 I =0» 85 

CH2 = XH2I*2/ (XH21*2+ ( 1-XH2I )*4 ) 

CHE=1 .-CH2 
CH=0* 

CH2P=0# 

CHP=0* 

CE = 0» 

READ (5*1) TINF*PREI *VEL0I *DENI *Z*CPI I 

1 FORMAT (6E10*3) 

DO 2 M= 1 *20 

DI (M)sDENI 
V I (M ) = VEL0I 
PINF (M)=PREI 

ETHP(M )=-l .AIZE-S/VI (M )**2 

CPL(M)=19.52E7 

QO (M )=0* 

2 CONTINUE 
KDD=1 

6 CALL SHOKLY 

IF (KDD#E0.2) GO TO 3 

IF ( (QOl ( 1 )-Q0( I ) )/Q01 ( 1 ) •LT.0.03) GO TO 3 
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QO ( 1 )=Q01 ( 1 ) 

QOl (MBP)=QOl (MA2 ) 

DO 5 Msl •MBP 
YI=0» 

S=M*1 • 

VELOIsABSlVELOl ) 

TSDIM=TS(M)*VEL0I**2/CPI I 

call PREC2 (M*Y1 ,TSDIM»PREI • ARFA(M)t QOl(M). VELO I *DEN I t TP*PPt VPt 
1 Ult VT« DP. CHM* CH2M. HE I .CPC) 

VI CM)=VT 
DI (M)=DP 
PINF (M)=PP 

ETHP(M)=HEI /VI(M)**2 
5 CONTINUE 
KDD=2 
GO TO 6 
3 CONTINUE 

WRITE (6.332) 

332 FORMAT < / • 3X .*TEMP* . 1 2X . *PRE* . 1 2X , *V- VELO* * 1 OX . *U- VELO* • 9X . *OEN* . 1 
! 2X.*CH*. 12X. *CH2+*. 12X.*ENTH*, 4X,*Y2*. 5X.*M*) 

DO 11 MB1.MA2 
Y1 =0. 

TSDIMsTS(M)#VELOI**2/CPI I 
9 VELOI=ABS(VELOI ) 

call PREC2 CM. Y1 .TSDIM.PREI . ARPA(M). QOl CM). VELO I . DENI . TP .PP. VP. 
1 UI. VT. DP. CHM. CH2M. ENTHC.CPC) 

Y2 = Y1/ (NS (M )#23 ) 

WRITE 16.331) TP. PP . VP. Ul . DP. CHM . CH2M . ENTHC . Y2.M 
331 FORMAT ( 2 X ♦ 8E 1 4 • 4 ♦ 3X . E 1 2 • 4 . 3X ♦ I 2 ) 

IF (Y1 .GT.10.0) GO TO 1 1 
IF (Yl.GT.3.0) GO TO 16 
IF (Y1 .GT.l .0 ) GO TO 15 
Y1 =Y1 +0.2 
GO TO 9 

15 Yl=Yl+0#5 
GO TO 9 

16 Y1=Y1+1 .0 
GO TO 9 

11 CONTINUE 
STOP 
END 


KD SUBROUTINE PREC2(S. Y. TS . PRE I . ZTA . QO, VELO I . DENI. TP. PP. VP. 
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I 


1 UI« VT*DP* CHM» CH2M» HE I » CPC) 

COMMON ARF(3) *SUM (12 )* I 1 ♦ T ( 50 ) • TT » Y I ( 3 ) * YD ( 3 ) « NA *ROH ( 3 ) * IERR( I 2 ) 
l«VF(3*2)f DH«CK *C«K *OQf Y1 .SI « E3 ( 3 ) . J . VE » DE 
DIMENSION V(50 ) .U(50) »P(50) .PHT(50) .PH (50) .CH2(50) .CH(50) .DEN (50) 
REAL NA . NB . NA 1 . MH2 . MHH . MH . MHE 
TT = TS 
QQ=Q0 
Y1 =Y 
SI =S 

VEL0I=-VEL01 

VE=VELOI 

DEsDENI 

ROH( 1 )s4.1E-18 
ROH(2 )=8*2E-18 
R0H(3)=2.1E-18 
DH=6«6256E-27 
CK=1 .38054E-16 
C=3.0E10 
VF( 1 . 1 )=8.7E15 
VF(2.2 )=3.75E15 
VF (3. 1 )=3.75E15 
VF(3.2 )=I .15E15 
VF(1 .2 )=5*02EI5 
VF(2. 1 )*5.02E15 
VI=VEL0I*SIN(ZTA) 

U I = VELO I ♦COS ( ZTA ) ♦ ( - 1 ) 

V(1 )=VI 
U( 1 )sUI 
P( 1 )=PREI 
DEN( 1 )*DENI 
CH2( 1 )*0# 

CH( 1 )=0. 

A314.8E12 

D*4.5El2/2. 

YI ( 1 )=1 .0 
YI (2 )*0*875 
YI <3 )»0* 

YD( 1 )»0«0 
YD(2)=0«125 
YD(3)*1 • 

R*8.3413E7 

PHT( 1 )= (VELO I *♦2 )/2*+l *527*R*145. 

PH( 1 )=1 .527^145.*R 
T< 1 )-145. 
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20 DO 1 U I»1 *49 
11 = 1 

NA=<7«2431 122E2240.85 )*P( 1 )/{T(l )*10*> 

NA1=NA»1 •OE-6 
ARF( 1 )=ROH( 1 >*NA1 
ARF(2 )=R0H(2 )*NA1 
ARF ( 3 ) =ROH ( 3 ) *NA 1 
DEN( I+l )=DEN1*VI/V( I ) 

V ( I + 1 ) = ( DEN I * ( V I **2 ) -P ( I ) +P ( 1 ) ) / ( DEN I * V I ) 

U( 1+1 )=U1 

CALL QRADIA (Y*S*QR) 

PHTI=1 •458*R*145.+(VEL01**2)/2. 

PHT(1+1 )=(0ENI*VI*PHT1“QR>/(DEN1*V1 ) 

PH ( I + l >=PHT( 1 + 1 )-(V( I + l )+*2 + U( 1+1 )**2 )/2* 

T ( I + l )= (PH( I+l )-(5./4**R*T ( 1 )+A/2. )*CH2< 1 ) - ( 3 • /4 • *R*T ( 1 ) + D)*CH( I ) ) 
l/( 1.458*R) 

CALL PCH2 <PCHI .PCHD) 

N0=l • 

CH( I+l )=NB*PCHD 
CH2( I+l )=NB*PCHI 

P < I + l ) = DEN( I + l )*R*T (1 + 1 )* ( I 76* 1 7/4 00* +0«5# (CH2 (I + l )+CH( I + l > ) ) 

IF (ABS( (V( I+l >“V( I ) )/V( I ) )*GT*0.01 ) GO TO 59 
IF (ABS( (T( I + l )-T( I ) )/T( I ) )*GT.0.01 ) GO TO 59 
IF (ABS( (P( I+l >-P( I ) )/P( I > >*GT.0*01 ) GO TO 59 
IF ( ABS ( (PH ( I + l )-PH ( I ) )/PH ( I ) ) .GT.O.Ol ) GO TO 59 
IF ( ABS ( (PHT( I+l )-PHT( I ) )/PHT( I ) ) •GT.0,01 ) GO TO 59 
IF (ABS( (CH( I + l )~CH( I ) )/CH( I ) ).GT*0*0l ) GO TO 59 
IF <ABS( (CH2( I + l )-CH2( I ) )/CH2( 1 ) )*GT«0,01 ) GO TO 59 
IF (ABS( (DEN( I+l )-DEN( I ) )/DEN( I ) ),6T«0,01 ) GO TO 59 
GO TO 75 
59 GO TO 111 
111 CONTINUE 

75 CHH=0.8019-CH2( 1 + 1 )-CH( I + l ) 

CHE=0, 1981 
CHM*CH( I+l ) 

CH2M=CH2 ( I+l ) 

HE I =PH( I+l ) 

TP=T( I+l ) 

PP=P( I+l ) 

VP»-V( I+l ) 

OP=DEN( I+l ) 

VT= (VP**2+UI**2 )**0*5 

RETURN 

END 
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SUBROUTINE QRADI A ( Y ,S « QR ) 

COMMON ARF ( 3 ) * SUM ( 12 ) « I *T (50 ) * TS* YI (3 ) • YD (3 ) *NA *ROH (3 ) * I ERR (12 ) * VF 
1 (3*2)*DH*CK*C<K1 .Q0*YY*SS«E3(3) 

COMMON /FFF/ Z 

D I MENS ION E2 ( 3 ) * ARY (3)»E1 (3)»V1 (3*2) 
external FX5*Fxa*FX3 
FI =6«256E-27/ ( 1 • 3805E- 1 6»TS ) 

VI (1 » 1 )=VF( 1 ♦ 1 )*F1 
VI (1*2) =VF ( 1 *2 )*F1 
VI (2* 1 )=VF(2* 1 )*F1 
VI (2*2 )=VF(2*2)*F1 
VI (3* I )=VF(3* 1 )*Fl 
VI (3*2)=VF(3*2)*F1 
EPS=0*002 
DO 50 K=1 *3 
K1 =K 

30 CALL ROMBS ( V 1 ( K , 2 ) * VI ( K * 1 ) * FX2 * EPS * SUM ( K+3 > * I ERR ( K+3 > ) 

IF ( IERR(K+3) .EQ.C) go to 50 

K3=K+3 

PRINT 31* K3* IERR(K+3) 

31 FORMAT (*K=** I2**IERR=** 13) 

50 CONTINUE 

QRl -0 

DO 100 K=l *3 

GR=0*5772 

ARY(K ) = ARF(K)*Y 

IF (Y*LT*0*001 ) GO TO So 

A0=0*26777343 

A1 =8«634 7608925 

A2=18*0590169730 

A3=8. 5733287401 

00=3*958469228 

81 =21*09965308 

82=25*63295614 

83=9*5733223454 

IF (ARY(K)*LT.1 . ) GO TO 200 

X = ARY (K ) 

El (K )=EXP (-X )* ( AO+A 1 *X+A2*X**2+A3*X»*3+X**4 )/ (X* ( BO+B 1 *X+B2*X**2+8 
13*X»*3 + X*»4 ) ) 

GO TO 220 
200 X=ARY(K) 

El (K ) =-0*57721 566-ALOG(X)+X-X**2*/4*+X**3*/18*-X»*4*/96*+X**5*/600 
1 * -X**6 * /4320 * + X **7* /35280 * “X»*8 * /322 560 * +X**9 « /3265920 * -X** 1 0 * /362 
288000* 



220 E2 <K )= (EXP<“X)-X*E1 (K ) ) 

E3(K )=(EXP(-X)-X*E2(K ) )/2. 

GO TO 35 
33 E2(K)=1* 

E3(K )=0.5 

35 QR4 = E3 ( K ) *SUM ( K+3 ) * 1 5 • *Q0 / ( 3 . 1 4 I 59**5 ) 

QR1= QR4+QR1 

100 CONTINUE 

0R=QR1*2.*3*14159 

RETURN 

END 

FUNCTION FX2(X) 

COMMON ARF(3) •SUMU2 ) » I . T ( 50 ) . TS • YI (3 ) . YD(3 ) * NA tROH ( 3 ) t lERR I 1 2 ) fVF 
1 (3*2) *DM*CK,C*K*QO 
DAT=5*6697E-5 
FX2=X**3/(EXP(X )-l • ) 

RETURN 

END 

SUBROUTINE PCH2 (PCHl.PCHD) 

COMMON ARF(3 ) *SUM (12)*I*T(50)*TS*YI (3)* YD (3) *NA*ROH (3) * IERR( 12) *VF 
1 (3f2)*DH*CK*CtK*90*Y*S*E3(3) *J1 *VEL0I *DENr 
DIMENSION VI (3*2) *CHl (3) «CH2(3) 
external FX4 

FI =6*256E-27/ ( 1 •3805E- 1 6*TS ) 

VI (1 * 1 >=VF( 1*1 )*F1 
VI (1 *2 )=VF( 1 *2)*F1 
VI (2* 1 )=VF(2*1 )*Fl 
VI (2*2)=VF(2*2)*F1 
VI (3* I )-VF<3* 1 )*F1 
VI (3*2)=VF(3*2)*F1 
EPS=0*001 
DO 59 J=1 *3 
JI =J 

CALL ROMBS ( V 1 ( J • 2 ) ♦ V I ( J * 1 ) * FX4* EPS* SUM ( J+9 ) * I ERR ( J+9 ) ) 

CHI (J)sYI ( J )*E3( J )*SUM ( J+9 )/ ( 1 .3805E-16*TS) 

CH2 ( J ) 3 YD ( J ) *E3 ( J ) *SUM ( J+9 ) / ( 1 • 3805E- 1 6*TS ) 

59 CONTINUE 

PDHI=CHl ( 1 )+CHl (2 )+CHl (3) 

PDHD=CH2 ( 1 )+CH2 (2 )+CH2 ( 3 ) 

XI -Q0*1 5**2**3*28E-24/ (DENI*VEL0I*3# 1 4159**4 )* (“1 • ) 

PCH1=PDHI»X1 
g PCHD=PDHD*X1 

S RETURN 

END 
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FUNCTION FX4(X) 

COMMON ARF (3 ) *SUM (12)*I*T(50)«TS*YI (3)* YD (3) *NA*ROH (3 ) • I ERR ( 1 2 ) • VF 
1 (3*2)«DH*CK*C*K*Q0,Y*StE3(3) «J 
FX4=X**2/(EXP(X )-l ) 
g RETURN 

END 


SUBROUTINE SHOKLY 

COMMON / BODY/ ZTA(20>« CK(20>* X(20)* DX * R ( 20 ) . BATA ( 20 ) * MUREF* EPS 
COMMON /CONST/ GAMA ♦ DY ( 5 1 ) * Y ( 5 1 > ♦ CC3 ( 20 . 5 1 ) *CCCl*ICONT 
COMMON /CONSTl/ AL • AK « CT * K . K 1 I ♦ MBP , MA 1 , M A2 

COMMON /C0NST3/ AH2 ( 3 • 7 ) t AHP ( 3 « 7 ) « AH2P ( 3 * 7 ) . AHE ( 3 * 7 ) . AHEP ( 3* 7 ) « 

I AE(3«7) *AH<3»7)*KE 

COMMON /C0NST6/ YH*YK*TB 

COMMON /DEPEND/ T(20*51 ).P(20»5I )tD(20.51 )*U(20»51 )»V(20#51 ) *XHE(5 

I I ) *XHH( 51 ) » XHI5I ) *XHP<51 ) *CP(51 ) *XHEP(51 ) »XHE2P(51 ) *XEMlN(5l ) 

1 «THP1 (51 ) •TH(51 ) 

COMMON /HEAT/ QCON f D I FU • RAT » QTOTL 
COMMON /INF/ VINFt D INF* CP I 

COMMON /INFI/ TINF* PREI* VELOI * DENI* M* ALT* 1111* CPII*XH21 
1 . norad *NOPRE 

COMMON /PERC/ V I ( 20 ) ,P I NF ( 20 ) * ETHP(20), O I ( 20 ) , CPL ( 20 ) * 001 (20)* 

I PINFN(20) 

COMMON /RAD/ COF 1(58*51)* QFXP ( 5 1 ) * QFXM ( 5 1 ) * QFX ( 5 1 ) • DQY ( 5 1 ) 

COMMON /Rl/ A(58)*B(58)*F1 (58)*F3(58)*F(5B)* VO ( 1 0 ) * W1 ( 58 ) 

COMMON /SDFDIS/ NS ( 20 ) ♦ NS 1 * NS2 ♦ PS 1 * AS ( 20 ) *NSS2 

COMMON /SHOCK/ DS (20 ) * TS(20)* VS(20)* US(20)* PS ( 20 ) * MUS ( 20 ) • TKS ( 
120 ) ,ARFA(20 ) *CPS(20 ) ,THS(20) 

COMMON /SPHT/ CH2 ( 5 1 ) * CH ( 5 1 ) « CHP ( 5 1 ) * CHE ( 5 1 ) 

COMMON /TRANS/ TK ( 5 1 ) * MU ( 5 1 ) * TKD I 

REAL MUREF* MUHH * MUHE * MAC I * MU* MUS *NS * NS 1 * NS2 *MUUl 

INTEGER HIST 

RAT=1 • 1 0 

MBP=1 1 

MAI sMBP+l 

MA2=M0P-1 

TB=4000 

CALL DATA IN 

IC0NT=1 

DO 43 13*1*20 

PINFN( I3)=PINF( I3)/(DI ( I3)*VI ( I3)**2 ) 

43 CONTINUE 
NS2=0# 
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K=51 
K1 1=K-1 
II I 1-1 
DIFU=0»005 
VELOl=ABS(VELOI ) 

CALL BODY 
10 M»1 

IF (1111 .GT*1 ) GO TO 2 

32 IF (M.EQ*1 ) GO TO 5 
NS(M )=NS(M-1 ) 

GO TO 15 

5 NS(1 )=0* I 

15 CALL SHOCK 
M = 1 

33 CONTINUE 
VINF=VI (M) 

DINF*D1 (M ) 

CPI=CPL (M ) 

IF (1111 *GT.l ) GO TO 81 
IF (M»EQ* 1 ) GO TO 6 
NS(M)sNS(M-l ) 

GO TO 16 

6 NS(1 )=0,1 

16 IF (1111 *GT. 1 ) GO TO 2 
ZTAA=ARFA(M) 

VENF=VINF*l .E-B 

UT*VENF#SIN(ZTAA )*(1 .+0.7476*( 1 .-XHEl ) ) 

CTU«-545 • 37+6 1 • 60B+UT-2 • 2459*UT **2+0 •039922*UT**3-’0 • 00035 1 48*UT **4 
I +0 *0000012361*07**5 

CHT=5 *661 -0*52661 *UT + 0 ,020376*07**2-0 *00037061 *UT**3+0* 0000034265* 
1 UT**4-0 * 0000000 1 2206*UT**5 
CH-CHT-0*3167*( I *-XH2l ) 

CT=CTU+61 *2*(1 *-XH2 I ) 

IF (M*GT* 1 ) GO TO 4 
DO 21 N=:l *51 
P(M*N) = 1 • 

N2=N-1 

U(M*N)=1 */50.*N2 

V(M,N)=1 */50**N2 

T(M,N)= (0*75/50* )*N2+0* 25 

TSSsT(M*N)*TS(M )*VINF**2/CPI 

PSS=P ( M * N ) *PS ( M ) *D I NF* V 1 NF**2 

TSSl = (PSS/1013250. )**AL 

DSS= 0.001 292* ( CT*TSS 1 /TSS )**(!* /AK ) 



106 


D ( M ♦ N ) = DS5/ < OS ( M ) #D I NF ) 

21 continue 

GO TO 8 

4 DO 22 N=1 *51 
T(M*N)=T(M-1 ,N) 

P(M,N)=P(M-1 *N) 

D(M*N)=D(M-1 *N) 

U(M,N)=U(M-1 *N) 

V(M*N)=V(M“1 *N) 

NS(M)=NS<M-1 ) 

22 CONTINUE 
GO TO 0 

2 CONTINUE 
CALL SHOCK 
M=1 

81 NST=NS(9) 

8 UST=U(M*20) 

ZTAA = ARFA (M ) 

UT = VENF*S1N(ZTAA )*(l •+0.7476*( 1 •-XH2I ) ) 

CTU= -545 .37+61 • 6O0*UT-c. . 2459*UT**2+ 0 • 039922*UT**3-0 • 00035 1 48*UT**4 
1 +0.0000012361*UT**5 

CHT=5 .661 -0.52661 *UT+0. 020376*UT**2 -0.0003706 1 *UT*»3+0. 0000034265* 
1UT**4“0. 0000000 12206*UT**5 
CH=CHT-0.3167* ( 1 .-XH2I ) 

CT=CTU+6l .2*( 1 .-XH2 I ) 

DST=D(M*20) 

TST=T(M*2 ) 

DEFU=DIFU 

HIST=1 

105 DO 23 Nl=l .51 
N=52-N1 

TSS=T(M.N)*TS(M )*VINF**2/CPI 
RHO=D< M.N )*DS<M )*D1 NF* 1000. 

PI 1=TSS*(RH0/<1000. *0.001292 > )**AK 
P12=1013250.*<P1 1/CT )**( 1 ./AL) 

P13=0. 1*P12 

call NOBDEN ( TSS « RHO ♦ H2 * H . HPLUS * HE * HEPLUS * HE2PLUS * EM I NUS * XH2 I ) 
AAl=l./< H2+H+HPLUS+HE+HEPLUS+ HE2PLUS+EM I NUS ) 

Rl =1 .98726 
EMINS-EMINUS 
ND= ICONT/2 
ND=ND*2 

IF (NORAD .EQ.l) GO TO 902 
IF (ND.NE. ICONT ) GO TO 902 
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IF (HIST.EQ.2) GO TO 902 

901 call ABSOCOF (TSS»RH0»P12*H2fH*HPLUS»HE*HEPLUS«HE2PLUS*EMlNS) 

903 DO 900 KCD=1*58 

COFl (KCD*N)=A(KC0) 

900 CONTINUE 

902 XHH(N )=H2*AA1 
XH(N)=H*AA1 
XHE(N)=HE*AA1 
XHP(N)=HPLUS*AA1 
XHEP (N ) *HEPLUS*AA I 
XHE2P ( N ) =HE2PLUS* AA 1 
XEMIN(N )=EMINUS*AA1 
Cl 1=R1*2.5 

IF (TSS*LT«6000* ) GO TO 51 

CPH2=R1 #(3.363+4.656E“4*TSS-5» 127E-8*TSS**2+2*802E-12*TSS**3-4*905 

1E“17»TSS**4 ) 

CPH = R1 * ( 2 .4751 64 +7* 366387E-5*TSS-2 • 537593E-8#TSS#*2+2 • 386674E-1 2*T 
lSS**3-4 .5514312-17*755**4 ) 

CEMIN=R1* (2.508-6 .332E-6*TSS+1 . 364E-9*TSS**2- 1 . 094E- 1 3*TSS**3+2 . 93 
14E-1B*TSS**4 ) 

GO TO 52 

51 CPH2=R1*(3. 100190+5. 1 1 1946E-4*TSS+5.26442E-8*TSS**2-3.490997E-1 1*T 
1SS**3+3.69453E-15*TSS**4 ) 

CPH=C1 1 
CEMIN=C1 1 

52 CPHE=C1 1 
CPHPL=C1 1 
CHEP=C1 1 

CX = CPH2*XHH (N)+CPHE*XHE< N)+CPH*XH (N >+CPHPL*XHP (N ) + CEM I N*XEM I N (N ) +C 
1HEP*XHEP(N) 

CX2=XHH ( N ) *2 • +XH ( N ) +XHE ( N ) *4 • +XHP { N ) +XEM I N ( N ) *0 • 5486E-3 
CX 1 =CX/CX2*4 . 1 81 E7 
IF (N.EQ.51 ) GO TO 9 
CP(N)=CX1/(CPI*CPS(M ) ) 

GO TO 23 
9 CPS(M )=CX1 /CPI 
CP(51 )= 1 . 

23 CONTINUE 

IF ( HIST.EQ.2) GO TO 905 
IF (NORAD .EQ.l) GO TO 905 
IF (ND.EQ. ICONT ) GO TO 904 
IF ( ICONT.NE.l . ) GO TO 905 
IF (M.NE. 1 ) GO TO 905 
DO 93 N4=l .51 
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QFX(N4 )=0* 

DQY(N4 )>=0» 

93 CONTINUE 
GO TO 905 

904 call HAOAT 

905 CALL TRANSP 

IF ( HIST*EQ#2) GO TO 83 
CALL ENERGY 
HIST=2 
GO TO 105 
83 CALL MOMENTM 

ICONTs ICONT+1 

IF (ABS( (U(M*20 )-UST )/UST)«GT. 0.040) GO TO 8 
IF < ABS ( (D(Mf20 )-DST )/DST) .GT.0.030 ) GO TO 8 
IF (ABS((T(M*2 )-TST )/TST) *GT. 0.015) GO TO 8 
IF (ABS ( (DIFU-DEFU)/DIFU) .GT.0.10) GO TO 8 
IF (ICONT .LT.4) GO TO 8 
WRITE <6»I13) M*NS(M )* ICONT *X(M) 

113 FORMAT (////•3X«*STATI0N M=*»I2*5X* *NS(M)=*» E12.4* 5X» *ITERAT10 
1 Ns*» I2*5X.*X(M )=**F4. 1 > 

WRITE (6*112) QCON.DIFU * QTOTL * ETHP ( M ) 

112 format (//*2X**QcON=**E1 1 .4*5X**DIFU = **E1 1 .4** QTOT = * • E 1 1 . 4 * 

1 * ETHP(M)=*,E1 1 .4 */ > 

WRITE <6.72) 

72 FORMAT ( / . 2X . *U-VELOC I TY* ♦ 6X . *V-VELOC I TY* . 8X . *PRE* . 1 1 X . *TEMP*» 1 OX 
l.*DENSITY#. 9X.* TO I M* . 9X ♦ *QFXP * . 6X . ♦QFXM* . 9X . *Y* * 1 OX * *N* , / ) 

DO 79 Nal ,51 
TT=T (M,N)*TS(M) 

TTl -TT#VINF»*2/CPI 

WRITE (6*75) U(M*N) *V(M.N) * P(M*N). T(M.N). D ( M *N ) * TT 1 .QFXP ( N ) • OFX 
1 M(N).Y(N).N 

75 FORMAT ( 1 X * 9E 1 4 .4 , 3X , I 2 ) 

79 CONTINUE 

WRITE (6.^2) 

62 FORMAT (//.8X. *XH2* . 1 2X ♦ *XHE* . 1 2X • *XH* , 1 2X . *XH+* . 1 1 X ♦ *XE-*. I 1 X . *C 
1P*.12X.* TK*. 12X,*MU*.5X.*N*,/ ) 

DO 65 N=1 .51 

TKKl =TK (N )*TKS (M)*CPI*MUREF 
MUUl =MU (N )*MUS (M)*MUREF 

WRITE (6.61) XHH(N). XHE ( N ) . XH ( N ) , XHP ( N ) . XEM I N (N > . CP (N ) , TKK 1 , MUUl 
I ,N 

61 FORMAT ( 1X.8E14.3.3X.I2) 

THPl (N )=TH(N) 

65 CONTINUE 



IF (NOPRE.EO.l) GO TO 24 
QOl (M)=QFXP(51 )*CCC1 
24 M=M+1 
123 1C0NT=1 

IF dill *£0.2 ) GO TO 66 
IF (M.LT.MAl ) GO TO 33 
GO TO 67 

66 IF (M.LT.MBP) GO TO 33 

67 ASS= (NS (4 )-NS( 1 ) ) /6 
AS( 1 )=NS( 1 ) 

AS(2 )=AS( 1 )+ASS 
AS(3 )-AS<2 )+2*ASS 
DO 7 M1=4,M8P 
90 AS(M1 >-NS(Ml ) 

7 CONTINUE 

IF dill .EQ* 1 ) GO TO 111 
GO TO 60 
111 NST=NS(3) 

1111=2 
GO TO 10 
60 CONTINUE 

IF (NORAD .EQ*1) GO TO 201 
IF (NOPRE.EQ.l) GO TO 201 
RETURN 
GO TO 202 

201 STOP 

202 END 
C 

C 

SUBROUTINE DATA IN 

COMMON /Rl/ A(58) ♦B<58)*F1 (58)»F3(58) •F(58)« VO ( 1 0 ) * W 1 ( 58 ) 

DATA (FI ( I )* 1 = 1 f58 )/ 14,595t 1 3,595. 1 3.21 5* 13.086» 13.016* 12.765* 

A 12.725*12.545* 12.284* 12.1 06*12.086* 12.082* 12.064*1 1 .884* 

B 1 1 . 196« 10.696* 10.246* 10.201*10. 1965* 10. 1955* 10.191* 10.146* 

C 9.696* 9.000* 7.000* 5.000* 4.000* 3.400* 3.020* 2.875* 

D 2.835* 2.655* 2.579* 2.552* 2.546* 2*519* 2.249* 1.988* 

E 1.898* 1.889* 1.887* 1.878* 1.788* 1.511* 1.131* .987* 

F .947* .850* .731. .668, .654* .591, .470* .396* 

G .315* .297* .216* 0.000/ 

DATA (F3( I ) * 1 = 1 *58)/ . 482E-3 * . 71 6E-3 * . 83 1 E-3 * • 879E-3 * . 900E-3 * 

A .934E-3*.966E-3*.992E-3*.105E-2*.l lOE-2, .1 1 3E-2 * . I 1 3E-2 * 

B . 1 14E-2* . 1 17E-2 * . 130E-2* . 153E-2* • 174E-2 , . 1 87E-2 * . 1 89E-2 * 

C . 1 89E-2 * . I 89E-2 * . 1 90E-2 * .205E-2 * .246E-2 * .403E-2 , .980E-2 * 

D .225E-1 * .400E-1 * .609E-1 * .782E-1 , .860E-1 * .gsgE”! * • 1 12E+0* 
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E • 1 18E+0* • 121E + 0* . 123E+0, • 149E + 0 + .273E +0, ,295E + 0» 

F «297E+0» .300E+0, *325E + 0 « *452E + 0f •905E + 0* • 1 70E+1 » .221 E + 1 » 

G .277E+1 * .410E+1 ,.587E+1 ♦.693E+1 « •833E+1 « . 138E+2* .250E+2» 

H .457E+2 * .699E+2 , • 125E + 3** 159E+4/ 

DATA (F ( I ) . 1 = 1 ,58 )/ 1 • 1 9 , 1 • 04 , • 986 , .967, .960, .948, .937, .929, .913, 

A .897, .890, .889, ,888, .881 , .849, .805, .770 , .752 , .750, .750, .750, 

B .748, .730 , .683 , .533 , .44 1 , .331 , .272, .2 36, .217, .210, .2 02, .192, 

C . 189, . 187, . 186, .1 75, . 156, . 143, • 1 39, • 139, • 138, • 135, . 121 , .09 72, 

D .0779, .0711, .0661 , .0581 , .0515, .0486, .0458, .0390, .0318, .0261 , 

E .0225, .01 89, .00794/ 

DATA (V0(I), 1=1,10)/ 10.196,12.084,12.745,13.051,1.888,2.549, 

A 2.855, .661 , .967, .306/ 

DATA (W1 ( I ) , 1=1 ,53)/ 35.4, 14.595, 13.595, 13.215,13.066,13.016,12.76 
A 5, 12.725, 12.545, 12.284 , 12.106,12.086, 12.032, 12.046, 1 1 .384, 

B 11.196, 10.696 , 10.246, 10.201,10.1 965, 10.1 955, 10.191,10.146, 

C 9. 696, 9 .000, 7. 000, 5. 000, 4. 000, 3. 4 00 ,3. 02 0,2. 875, 2. 835, 2. 655, 

D 2.579,2.552,2.546,2.519,2.249,1 .983,1 ,898,1.889,1.887,1 .878, 

E 1 . 788, 1.511,1.131 ,0.987,0.947,0.850,0.731 ,0.688,0.654 ,0.591 , 

F 0.470,0.396,0.315,0.297,0.216/ 

RETURN 

END 

SUBROUTINE ENTHAlP ( TEMP, CH2 , CH , CHE , CH2P , CE , CHP , CPC , ENTHC ) 
COMMON /ENTH/ THH2, THH, THHE, THHP, THE 

COMMON /CONST3/ AH2 ( 3 , 7 ) , AHP ( 3 , 7 ) , AH2P ( 3 , 7 ) , AHE ( 3 , 7 ) , AHEP ( 3 , 7 ) , 

1 AE(3,7) ,AH(3,7),KE 

IF (KE.NE. 1 ) GO TO 15 

DATA (AH (1,1 ),!=! ,7)/2.5,0.,0.,0.,0. , 2 . 547 1 62E4 , -4 . 60 1 176E-1/ 

DATA (AH(2,I),I=1,7)/2.5,0,,0.,0.,0, , 2 . 547 1 62E4 , -4 . 60 1 176E-1 / 
data ( AH(3, I ) , I = 1 ,7 )/2. 4751 64, 7. 366387E- 5, -2.5375936-8,2. 386674E-1 
1 2,-4.55l43lE-l 7 , 2 . 52 3626E4 , -3 . 749 1 37E~ 1 / 

DATA ( AH2( 1 , 1 > , 1=1 ,7 )/3. 057445, 2. 676520E-3, -5. 80991 6E-6, 5.521 039E- 
1 9,-1 .8l2273E-12,-9.889047E2,-2.299705/ 

data ( AH2 (2, I ),I=1,7)/3.100190,5.111 946E-4 , 5 . 26442 1 E-8 , -3 . 490997E- 
1 1 1 ,3.694534E-15,-B.773804E2,-1 .962942/ 

data (AH2(3, 1 ) , 1=1 ,7 )/3. 3630 0,4 .656006-4,-5. 1 2700E-8 , 2 • 80200E- 1 2 , - 
I 4. 905E-1 7,-1 .0180063,-3.716/ 

data ( AHP( 1 ,I ), I = l,6)/2.5,0.,0.,0.,0.,1 . 840334 E5/ 

DATA (AH2P( 1 , I ) , I =l ,6 )/2. 8 17375, 3. 6576 lOE-3, -7. 9655486-6,8.26 14 OE- 
1 9, -3. 090228E-12, 1 .8002738965/ 

DATA (AH2P(2, I ) , 1=1 ,6)/ 3.32817562,2.505067816-4, ) . 42245209E-7', -4 . 
1 45902420E-1 1 , 3 . 733756E- 1 5 , 1 .799747065/ 

DATA (AHE (1 ,I ), I=l,7)/2.5,0.,0.,0.,0. , -7 . 453749E2 , 9 . 1 53488E- 1 / 

DATA ( AHEP( 1 , I ) , 1 =l , 7 ) /2 . 5 , 0 .0 , 0 . 0 , 0 . 0 , 0 . 0 , 2 . 853426E5 , 1 .608404/ 
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DATA (AE (1*1), 1=1*7) /2«5*0*,0.»0.,0, , ~7 • 453750E2 * - 1 .173402E1/ 

DATA (AE(3* n * 1 = 1 .7 ) /2 . 508 * -6 • 3320E-6 * 1 .3640E-9*-! *094E- 1 3 *2 *9340E 
1 -18*-7.450E2*-1 •208E1/ 

DO 1 1=1*6 
AHP(2* I )=AHP( 1 * I ) 

AHP(3* I )=AHP( 1 * I ) 

AH2P(3* I )=AH2P(2, I ) 

AHE(2* I )=AHE( 1*1) 

AHE(3* I )=AHE( 1 * I ) 

AHEP (2*1) =AHEP (1*1) 

AHEP(3* I )=AHEP (1*1) 

AE (2* I )=AE( 1,1) 

1 CONT I NUE 
15 R=1 ,98726 

IF (TEMP. GT, 1000, ) GO TO 2 
J=1 

GO TO 10 

2 IF (TEMP.GT«6000, )G0 TO 3 
J = 2 

GO TO 10 

3 J=3 

1 0 CPH=R* ( AH ( J * 1 )+AH ( J»2 )*TEMP+AH ( J*3 )*TEMP**2+AH ( J *4 ) *TEMP**3+AH ( J * 5 
1 )*TEMP**4) 

CPH2=R<fr ( AH2( J* 1 ) + AH2 t J*2 )*TEMP+AH2( J*3)*TEMP**2+AH2( J*4 )*TEMP**3+ 

1 AH2( J*5 )*TEMP*#4 ) 

CPHP=R» ( AHP ( J* I )+AHP ( J,2 )*TEMP+AHP( J *3 )*TEMP**2+AHP( J*4 )*TEMP#*3+ 

1 AHP( J*5 )*TEMP**4 ) 

CPH2P=R* ( AH2P ( J * 1 )+AH2P ( J * 2 ) *T£MP+AH2P ( J * 3 ) *TEMP*»2+AH2P ( J * 4 ) *TEMP 
1 **3+AH2P ( J,5 )*TEMP**4 ) 

CPHE = R* ( AHE ( j* I )+AHE ( J*2 ) *TEMP+AHE ( J * 3 ) *TEMP**2+AHE ( J*4 )*TEMP**3+ 

I AHE ( J*5 )*TEMP**4 ) 

CPE=R*( AE( J* 1 )+AE ( J*2 )*TEMP+AE ( J*3 )*TEMP**2+AE ( J *4 ) *TEMP**3+AE ( J * 5 
1 )*TEMP*#4) 

R1 =R*TEMP 

THH=R1*(AH( J* 1 )+AH( J*2 )/2,*TEMP+AH( J,3)/3.*TEMP**2+AH( J*4 )/4.*TEMP 
1 **3+AH ( J*5 )/5»*TEMP**4,+AH( J*6 )/TEMP ) 

THH2 = R1 * ( AH2 ( J * 1 )+AH2( J*2)/2,*TEMP+AH2( J*3)/3.*TEMP**2+AH2(J,4)/4, 
1 *TEMP##3+AH2 ( J*5 )/5**TEMP**4+AH2 ( J*6 )/TEMP ) 

THHP=R I * ( AHP ( j * 1 ) +AHP ( J * 2 ) /2 . *T£MP+AHP ( J * 3 ) / 3 . ♦TEMP**2+AHP ( J * 4 ) /4 • 
1 *TEMP**3+AHP ( J*5 ) /5,*TEMP**4+AHP ( J * 6 ) /TEMP ) 

THH2P=R1*( AH2P( J* 1 )+AH2P( J*2 ) /2 , *TEMP+AH2P ( J * 3 ) /3 , «-TEMP**2+AH2P ( J « 
1 4 )/4.*TEMP**3 + AH2P( J*5 )/5«*TEMP**4+AH2P ( J*6 )/TEMP ) 

THHE = R1* ( AHE ( J* 1 ) +AHE ( J * 2 ) /2 ,*TEMP+AHE ( J * 3 ) /3 . *TEMP**2+AHE ( J*4 )/4, 
1 *TEMP**3+AHE( J*5 )/5,*TEMP**4+AHE ( J*6 )/TEMP) 
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THE=R1 * ( AE< J * 1 )+AE ( J *2 > /2 • *TEMP+AE ( J * 3 ) /3 • *TEMP**2 + AE { J . 4 )/4*TEMP* 
1 *3 -»-AE(J*5)/5.*TEMP**4 + AE( J*6)/TEMP) 

CPsCPH*CH+CPH2*CH2/2 .+CPHP*CHP+CPH2P*CH2P/2.+CPHE*CHE/4 .+CPE*CE/0. 
1 5486E-3 

CPC=CP/<2.3901E-8 ) 

ENTh=THH*CH+THH2*CH2/2*+THHP*CHP+THH2P*CH2P/2.+THHE*CHE/4 •+THE*CE/ 
1 0*5486E~3 

ENTHC=ENTH/(2.3901E-e ) 

C 

C CP**CAL/GM K CPC* .CM2/SEC2 K 

C ENTH**CAL/GM ENTHC* *CM2/SEC2 K 

C CH2* *CH* *CH2P* MASS FRACTION 

RETURN 

END 


SUBROUTINE BODY 

COMMON / BODY/ ZTA(20)* CK(20)* X(20)* DX » R ( 20 ) * BATA ( 20 ) * MUREF » EPS 
COMMON /CONST/ GAMA » DY ( 5 1 ) » Y ( 51 ) » CC3 ( 20 * 5 I ) 

COMMON /INFI/ TINF* PREI» VELOi . DENI* M. ALT* Mil* CPII*XH2I 

COMMON /INF/ VINF. D INF, CP I 

COMMON /HEAT/ QCON * D I FU * RAT * QTOTL 

REAL MUREF* MUHH , MUHE * MAC I * MU * MUS * MUH 

WRITE (6*11) 

11 FORMAT (/,5X**ZTA**9X*»B0DY CURVE** 6X * *W* * 9X * *MUREF* * 7X * *EPS* * 

1 8X**M»*/) 

RN=23* 

TREF=VEL0I**2/CPI I 

MUREF = 0.66E-6* (TREF**1 . l>/ ( TREF + 70 #5 ) )»10, 

EPSsMUREF/ ( DEN I *VELO I *RN ) 

EPS=EPS**0.5 
PI=3. 14159 
X ( 1 )=0* 

ZTA ( 1 )=PI/2* 

BATA (1 )=0* 

CK ( 1 ) = 1 . 

R( 1 )=0* 

DX = 0* 1 
A=l • 

Ml =20 

DO 5 M=2*M1 

7 FI =DX*SQRT( (A**2+R(M-1 )**2 )/ ( A**2+2 • *« ( M~ 1 )**2 ) ) 

F2 = DX*SQRT( (A**2+ (R(M-1 )+Fl/2. ) **2 ) / ( A**2+2 • * ( R ( M- 1 )+Fl/2* )**2) ) 

F3 = DX*SQRT( (A**2+(R(M-1 )+F2/2. ) **2 ) / < A**2+2 •* ( R ( M- 1 )+F2/2. )**2) ) 
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F4=DX*SQRT( (A**2+(R(M-1 )+F3)**2 ) / ( A**2+2 • * ( R ( M-1 >+F3)**2) ) 

R (M)=R( M-l )+ (FI +2«*F2+2.*F3+F4 )/6. 

W=SQRT( A**2+R(M)**2 ) 

CK(M)=A**2/(R(M)**2+W**2)**1 .5 
ZTA(M)=ATAN(W/SQRT(W**2-A**2 ) ) 

BATACM )=PI/2*-ZTA (M ) 

X<M)“X(M-1 )+DX 

WRITE (6*89 >ZTA (M ) » CK ( M ) ♦« ( M ) *MUR£F «EPS* M 
89 FORMAT ( 3X* 5E1 2 *4 * 4X * I 2 ) 

5 continue 

YY=l .00*( 1-RAT >/( 1-RAT**50 ) 

DY ( 1 >=YY 
Y( 1 )s=0. 

DO 99 N=1 *49 
OY(N+l )*RAT#DY(N) 

Y(N+1 )sY(N)+DY(N) 

99 CONTINUE 

Y(51 )sY (50 )+DY (50 ) 

DX-0* 1 
RETURN 
END 
C 
C 

SUBROUTINE SHOCK 

COMMON / BODY/ ZTA(20)* CK(20)* X(20>* DX * R ( 20 ) ^BATA (20 ) * MUREF* EPS 
COMMON /CONST I / AL* AK * CT * K *K i 1 * MSP • MA I 

COMMON /C0NST3/ AH2 ( 3 * 7 ) ♦ AHP ( 3 , 7 ) * AH2P ( 3 * 7 ) * AHE ( 3 ♦ 7 ) * AHEP ( 3 * 7 ) * 
1 AE(3.7) *AH(3*7)*KE 

COMMON /INFI/ TINF* PRE I * VELOI * DENI* M* ALT* 111 I* CPII*XH2I 
COMMON /PERC/ VI (20 ) *PINF (20 ) * ETHP(20). DI (20 ) * CPL ( 20 ) * Q0i(20)* 

1 PINFN(20) 

COMMON /SDFOIS/ NS ( 20 ) * NS 1 * NS2 ♦ PS I . AS ( 20 ) 

COMMON /SHOCK/ DS(20)* TS(20)* VS(20)* US(20)» PS ( 20 ) , MUS ( 20 ) * TK3 ( 
120 ) *ARFA(20 > *CPS(20 ) *THS(20) 

REAL MUREF* MUHH * MUHE * MAC I * MU* MUS *NS *NS 1 * NS2 
DIMENSION HS(20 )*HSS(20) ,EH(20) * DNS (20) 

HX2I= XH2I 

MBB=MBP-1 

W=8.3143E7 

AL=0*67389-0.04637*ALOG(XH21 ) 

AK=0»65206-0*04407*ALOG(XH2I ) 

AM=0*95252-0.1447*ALOG(XH2I ) 

AN=0.97556-0*16149*AL0G(XH2I ) 

WRITE (6*39) 
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39 format <9X**PS*f 15X»*DEN*» 15X**TEMP** 14X«*HS** 1 6X * *US* « 1 6X ♦ *VS* * 

1 10Xt*M**8X**T (K)*,/) 

IF ( I 1 I 1 •GT. I ) GO TO 1 
DO 1 1 I =l *20 
ARFA ( I ) =ZTA ( I ) 

11 CONTINUE 
GO TO 4 

1 DO 12 1=2* MSB 

IF (I .EQ-MBP ) GO TO 2 

DNS( I )s (AS( I + l )-AS( I-l ) )/(2.*DX) 

GO TO 30 

2 DNS(MBP)= (AS<MBP)-AS(MBB) )/DX 

30 ARFA(I)= ZTA( I )+ATAN{DN5{ I )/( 1 *+CK< I )*NS( I ) ) > 

12 CONTINUE 
ARFA( 1 )=ZTA( 1 ) 

4 CONTINUE 
VINF-VI (M ) 

DINF =DI (M ) 

CPI=CPL (M ) 

VENF=V I NF*1 .E-5 
ZTAA=ARFA (M ) 

UT=VENF*SIN(ZTAA )* ( 1 *+0*7476»( 1 *-HX2I ) ) 

CTU=-545* 37+61 . 608*UT-2 • 2459*UT**2+0 • 039922*UT**3-0 • 00035 1 48»UT**4 
1 +0*0000012361 #UT*-»5 

CHT=5* 661 -0,52661 *UT+0 , 020 376*UT**2-0 • 0003786 1 *UT* *3+0 • 00 00034265* 
1 UT**4-0, 0000000 1 2206*UT**5 
CH=CHT-0,3167* ( 1 ,-HX2l ) 

CT = CTU+61 ,2* ( 1 *-HX2 I ) 

IF (M,EQ* 1 ) GO TO 3 
DS1=DS(M-1 ) 

DS2=DS1 +1 • 

27 D9=DS2 

EHl =ETHP<M ) +SIN ( ARFA (M ) )**2* ( 1 *-l ,/DSl**2 )/2* 

P1=PINFN(M> +SIN(ARFA(M) )**2*( 1*-1*/DS1 ) 

GO TO 301 

9 EH2=ETHP(M) +S I N ( ARFA ( M ) ) **2* ( 1 • - 1 • /DS2**2 ) /2 • 

P2=PINFN(M ) +SIN( ARFA (M ) )**2*( 1 .-1 ,/DS2 ) 

GO TO 302 

3 DS1~6 
DS2 = 9 

28 D9=DS2 
NS1=NS( 1 ) 

EHIsETHPIM) +0,5*< 1 ,-l */DS1**2 ) 

PP=PINFN<M) +( I .-1 ,/DSl ) 
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PP1=BATA(M)**2*( ( 1 .*1 •/OSl )*< 1 •-NS2/( 1 .+NS1 ) )**2 > 

PI =pp-ppi 

301 OSSsDSl*DINF 
PSSspl »DINF*VINF**2 
XA=PSS/10 13250. 

YA=DSS/0.001292 
TSS=CT*- (XA**AL/YA**AK ) 

RH0=DSS*1 000 

CALL N08DEN (TSS « RHO. H2 » H .HPLUS . HE « HE PL US ♦ HE2PLUS.EM INUS. 

1 XH2I ) 

CALL VOLMAS (H2.H* HPLUS.HE. HEPLUS . HE2PLUS » EM I NUS « YH2 » YH. YHP ♦ 
1 YHE* YHEP.YHE2P. YE) 

YH2P=0 

KE=1 

CALL ENTHALP (TSS ♦ YH2 » YH. YHE* YH2P. YE« YHP. CPC ♦ ENTHC ) 

KE=2 

EHHl =ENTHC 

DH1=EHH1/ (VINF#*2 )-EHl 
IF (M.NE. 1 > GO TO 9 

19 EH2=ETHP(M) +0.5# U .-Z ./DS2#*2 ) 

PS1=PINFN(M) +(l.-l./DS2) 

PP1=BATA(M)**2* ( ( 1 .-1 ./DS2 )* ( 1 .-NS2/ ( 1 .+NS1 ) >**2 ) 

P2sPSl -PPl 

302 PSS=P2»DINF*VINF#*2 
DSS=DS2*DINF 
XAspsS/1013250. 

YA=DSS/0.001292 
Z=W#273. 15/2 .304 
TSSsCT* (XA**AL/YA**AK ) 

RHO=DSS*l 000 

call NOBOEN (TSS « RHOf H2 . H .HPlUS . HE . HEPLUS. HE2PLUS * EM I NUS . 

1 XH2I ) 

call VOLMAS (H2.H. HPLUSiHE. HEPLUS . HE2PLUS . EM I NUS » YH2 . YH. YHP. 
1 YHE. YHEP.YHE2P. YE) 

call ENTHALP (TSS . YH2 . YH . YHE ♦ YH2P . YE . YHP . CPC . ENT HC ) 

EHH2=ENTHC 
EH3=EHH2/ (VINF#»2 ) 

IF ( ABS( (EH3-EH2 )/EH2) .LT. 

DH2=EH3-EH2 

SL0P=(DH2~0H1 )/ (DS2-DS1 ) 

D3=DS2-DH2/SLOP 
OHl =DH2 

IF (D3.LT.5. ) GO TO 8 
IF (D3.LT.14.) GO TO 10 


0.005) GO TO 100 
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8 OS2=D9+O*l 

IF (M.NE* 1 ) GO TO 27 
GO TO 28 
18 DS2=D3 

IF (N.NE. 1 ) GO TO 9 
GO TO 19 
100 PS(M)=P2 
DS(M)=DS2 

TSS=CT* <XA*»AL/YA**AK ) 

TS(M > = TSS*CPI/ (VINF**2 ) 

EH(M >=EHH2/ < VINF*VINF ) 

VSP=~SIN<ARFA(M > )/DS(M ) 

USP=COS ( ARFA (M ) ) 

US(M )=USP*SIN(ARFA (M )+BATA (M ) ) +VSP*COS ( ARFA ( M ) +dAT A ( M ) ) 
VS(M)=“USP*C0S(ARFA(M)+BATA(M) )+VSP*SIN( ARfA (M )+BATA (M > ) 

THS(M ) = EH(M )+US(M )**2/2 

WRITE (6<29) PS(M)* dS(M)« TS(M)* EHCM), USCM)^VS(M)» M*TSS 
29 FORMAT ( 4X*6E17.4«4X* I2»4XiEl 1 .4 ) 

M = M+1 

IF (M.GT. MBP) GO TO 69 
GO TO 4 
69 RETURN 
END 
C 

SUBROUTINE NOBDEN (TSS*RHO* H2*HtHPLUS* HE » HEPLUS f HE2PLUS « EM I NUS t 
1 XH2I ) 

XH2=XH?I 
XHEsl .-XHa 
ALPHAH=2«*XH2 
ALPHAHE-XHE 
ALPHAC=0. 

T = TSS 

AM0=1 •008*ALPHAH+4.003*ALPHAHE+12*01 1 *ALPHAC 
R=8,31441E3 

compute EQUILIBRIUM CONSTANTS *- NATURAL LOGS 
ANALOG ( T ) 

STARN=6«02252E26* (RHO/AMO ) 

0*1 ./T 

C 1 =69 • 939357-5 1 964 . »B 
C2=-49* 234384- ( 1 • 5*A ) + 1 578 1 O • *8 
AKK=ALPHAHE»STARN 
IF( ALPHAHE.EQ.O* )G0T04 
C3=-50»620676- ( 1 • 5*A ) +285287 . *B 
C4=-99« 161915- (3.*A )+9l6687**B 
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IF<C4#GT#741 • )C4 = 741 • 

IF(ALPHAC~0. )5,7»5 

05=292.632558- ( 1 .25*A )- 1 97547.4*6 
C6=2l 8. 091 51 2-A- 146258.4*6 
C7=144.31 1945-( .5*A ) -102732.0*6 
08=68.629830-40279 .5*6 
C9=-49. 522066-1 .5*A+ 1 30774 .8*8 
01 0=-97. 657838-3. *A+4 13782 .3*6 
FORMAT ( IH //) 

OOMPUTE SPEOIES NUMBER DENSITIES 
AK=ALPHAH*STARN 
AK1=EXP(C1 )/(4.*AK) 

AK2=2.*AK*EXP(C2) 

El =1 .-AKl* <SQRT ( 1 . + (2. /AKI ) )-l . ) 

E2=( 1 ./AK2 )*(SQRTU .+2.*AK2)-1 . ) 
IF(ALPHAHE.EO.O. )G0T08 
AK3=AKK*EXP(03) 

AK4=AKK*EXP (04 ) 

AK34=EXP(C3)/AK4 
A1 =ALPHAH/ALPHAHE 
E2AL=E2*A1 
D= ( I ./AK3 )+ (E2AL ) 

A2=l .+E2AL+AK34 

E4=.5*A2* (SORT ( I .+( (4.*AK34 )/(A2**2 )))-!•) 

E3=.5*D*(SQRT( 1 .+ < ( 4 ./AK3 ) /D**2 ) )-l • )-E4 

E12=E3+2.*E4 

I F ( ALPHAC-0 • ) 9 ♦ 99 . 9 

AKA=AK*SQRT (2.*E1 *AK1 ) 

AK5= ( AKA**4 )/EXP ( 05 ) 

AK6= (AKA**3 )/EXP(C6 ) 

AK7= ( AK A**2 )/EXP ( 07 ) 

AK8=(AKA)/EXP(08) 

acnstr=alphao*starn 
AK9= ACNSTR*EXP ( 09 ) 

AKl 9= I ./AK9 

EA=1 ./ ( 1 .+AK5+AK6+AK7+AK8 ) 

E5=AK5*EA 

E6=AK6*EA 

E7=AK7*EA 

E8=AK8*EA 

AA=(4.*E5 )+ (3.*E6 )+ (2.*E7)+E8 
IF(ALPHAHE.EQ.O. )G0T01 0 

El 1 = { (E2*ALPHAH )+ ( (E12 )*ALPHAHE) )/ALPHA0 
E10=l •/( 1 .+( (E2*AK)+ (E12*AKK ) ) *EXP C 0 1 0-09 ) ) 
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E9=*5*(E1 1+AK19 >♦ (SQRT( 1 •+ ( ( 4 • *AK 1 9 ) / ( E 1 1+A< 19 )**2 ) )-l . )“E1 0 
G0T099 

10 E10=l./(1 •+(E2*AK )*(EXP(C10-C9) ) ) 

EE=E2*< ALPHAH/ALPHAC )+AK19 

E9=.5*EE* (SORT ( 1 . +4 • *AK 1 9/ ( EE*EE ) )-l • )-ElO 
99 CONTINUE 

H2=*5*E1*AK 

H= ( 1 .-E1-E2 )*AK 

HPLUS=E2*AK 

IF(HPLUS.LT.0* )HPLUS=0, 

IF ( ALPHAHE*EQ*0. )G0T01 1 
HE= ( 1 .-E3-E4 >*AKK 
HEPLUS=E3*AKK 
HE2PLUS-E4*AKK 
EMINUS=HPLUS+ (E12*AKK ) 

Z=((l*-(£l/2. )+E2 )*AlPHAH) + ( ( 1 *+El2 )*ALPHAHE ) 

G0T050 

11 Z= < 1 •“(El/2« )+E2 >*ALPHAri 
HE = 0* 

EMINUS=HPLUS 

HEPLUS-0* 

HE2PLUS=0. 

50 IF(ALPHAC-0, )12* 14, 12 

12 H2=H2-< (El/( 1 .+E1 ) )*AA*ACNSTR) 

H = H- ( < { < 1 .-El )/ < 1 . + E1 ) )*AA )- ( E9* < ( 1 .-E2 )/ (2*~E2 ) ) ) )*ACNSTR 
HPLUSaHPLUS- (E9*( ( 1 .-E2 )/ (2.-E2 ) >*ACNSTR ) 

IFIHPLUS.LT.O. )MPLUS=0. 

IFCALPHAHE.EQ.O. )G0T077 

BBB= (E3*( 1 .-E3 )* ( 1 • +E 1 0 ) ) / ( (E3* (2.-E3 ) ) + Ai ) 

CC=(2#*E4*( 1 .-E4 ))/((! . + (2#*E4 )-(E4**2 ) )+Al ) 

HE=HE+ ( BBB»ACNSTR ) 

HEPLUS=HEPLUS- ( ( BBB-CC ) *ACNSTR ) 

HE2PLUS=HE2PLUS- ( CC* ACNSTR ) 

EMINUS=EMINUS+ { <E9/ (2.-E2 ) )+ (2**E10 )-BBB-CC )* ACNSTR 

Z = Z+ ( 1 .-AA/ ( 1 .+E1 )+E9/( 2.-E2 ) +2 . *E 1 0- ( E3* ( 1 .-E3 )* ( 1 .+E10 ) )/ (E3 + A1 ) 
1-CC)*ALPHAC 
G0T013 
77 HE=0. 

HEPLUS=0. 

HE2PLUS=0. 

D0=E9/(2.-E2 )+2**ElO 

EMINUS=EMINUS+DD*ACNSTR 

Z=Z+ ( I .-AA/ ( 1 .+E1 )+DD)*ALPHAC 

13 CH4=E5*ACNSTR 



119 


CH3=E6*ACNSTR 

CH2=E7»ACNSTR 

CHsE8»ACNSTR 

CN= ( 1 .-E5-E6-E7-E8-E9-E1 0 ) *ACNSTR 
CPLUS=E9*ACNSTR 
C2PLUS=E1 0*ACNSTR 
G0T088 
14 CH4=0. 

CH3sO* 

CH2sO« 

CH = 0* 

CN=0* 

CPLUS-0* 

C2PLUS=0# 

88 CONTINUE 

39 RETURN 
END 
C 

subroutine TRANSP 

COMMON / BODY/ ZTA(20)f CKC20)* X(20)* DX ♦ R ( 20 > « B ATA (20 ) * MUREF « EPS 
COMMON /CONST/ GAMA * DY ( 5 1 ) f Y ( 5 1 ) » CC3 ( 20 * 5 1 ) 

COMMON /CONST I / AL » AK ♦ CT 

COMMON /DEPEND/ T(20,51)*P(20,51 ),D(20«51 )*U( 20*51 )*V(20«5l )*XHE(5 
1 1 ) *XHH(51 ) * XH(51 ) .XHPC51 ) *CP(5I ) .XHEP(5l ) *XHE2P(5I ) *XEMIN(51 ) 
COMMON /INF/ VINF* D INF* CP I 

COMMON /INFI/ TINF* PWEI* VELOI* DENI* M* ALT* 1111* CPI I 
COMMON /SDFOIS/ NS ( 20 ) * NS 1 * NS2 

COMMON /SHOCK/ DS(20)* TS(20). VS(20)* US(20)* PS ( 20 ) , MUS ( 20 > * TKS ( 

120 ) ,ARFA(20 ) 

COMMON /TRANS/ TK <51 ) * MU ( 51) * TKD I 

REAL MUHH* MUH* MUHE * MUD* MUREF * MSD* MU* MUS 

K»51 

DO 10 Nl=l*51 

N=52~N1 

NNN=N 

TEMP-T(M*N)*TS(M)*VINF**2/CPI 

MUHH = 0*66E-6*TEMP*»1 • 5/ < TEMP+70. 5 ) * 1 0* 

MUHE=1 .SSE-S* <TEMP)**1 • 5/ ( TEMP+97.8 ) * 1 0 . 

TKHH = 3.21 lE-5+5«344E-8*TEMP/ (6*718E-2 ) 

TKH s=2*496E“5+5* 129E-8*TEMP/ <6#718E-2 ) 

MUH= TKH*4./(15»*1 .987) 

TKHE=MUHE*15.*1 •987/16* 

PCI 1-1 • 

PC22=1 • 
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PC33=1 • 

PC12=( 1 •+(MUHH/MUH)**0«5*(0«5)**0*25)**2/SQRT(24. ) 

PC13*( 1 •+SQRT(MUHH/MUHE)*2***0*25)*»2/SQRT( 12* ) 

PC21=( 1 .+SQRT (MUH/MUHH )*2.**0,25 )**2/SQRT ( 12* ) 

I- PC23=(1 *+SQRT (MUH/MUHE )*4«**0«25 )**2/SQRT ( lO. ) 

o PC31 = ( I .+SQRT (MUHE/MUHH)* (U.tjU ) **U .ciij ) **^;/bQR l ) 

PC32= < 1 •+SQRT (MUHE/MUH )*(0*25)**0.25 )**2/SQRT (40. ) 

PCHH=XHH < NNN ) *PC 1 1 +XH ( NNN ) *PC 1 2 + XHE ( NNN ) *PC 1 3 
PCHsXHH ( NNN ) *PC2 J +XH ( NNN ) *PC22 + XHE ( NNN ) *PC23 
PCHE*XHH { NNN ) *PC31 +XH ( NNN ) *PC32+XHE ( NNN ) *PC33 
PC=PCHH+PCH+PCHE 

MUD-(XHH(NNN)*MUHH)/PCHH + ( XH ( NNN ) *MUH ) /PCH + ( XHE ( NNN ) *MUHE ) /PCHE 

TKDs (XHH (NNN)*TKHH ) /PCHH + ( XH ( NNN ) *TKH ) /PCH + ( XHE ( NNN )*TKHE )/PCHE 

TKDsTKD*4*18E7 

IF (N.GT. I > GO TO 5 

TKDI =TKD 

5 IF (N.NE.K) GO TO 2 
MSD-MUD/MUREF 
TKSD=TKD/(NIUREF*CPI ) 

MUS(M ) sMSD 
TKS<M ) sTKSO 
2 CONTINUE 

MU (N ) -MUD/ ( MUREF*MSD ) 

TK (N ) =TKD/ ( MUREF*CP I *TKSD ) 

10 CONTINUE 
RETURN 
END 
C 
C 

subroutine ABSOCOF (TSStRH0.Pl2*H2»H* MPLUS*HE *HEPLUS* HE2PLUS* E 
IMINS ) 

*** CALCULATION OF ABSORPTION COEFFICIENTS *** 

COMMON /Rl/ A ( 58 ) *B < 58 ) «F1 < 58 ) f F3 (58 ) »F ( 58 ) » VO(IO) 

DIMENSION S(10)*G(10) 

REAL NH2 » NH f NHP *NHN * NE « NHE * NHEP * NHE2P 
NH2=H2/ 1 000000 • 

NH=H/1 000000. 

NHP=HPLUS/ 1 000000 • 

NHNaO. 

NHEssHE/ 1000000. 

NHEP=HEPLUS/1 000000. 

NHE2P=HE2PLUS/ 1 000000 . 

NE=EM I NS/ 1 OOOOOO . 
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P=P12 

T=*TSS 

ST s I57780*/T 

■»** FREE-FREE AND BOUND-FREE TRANSITIONS FOR HYDROGEN *** 

Cl * 9#93E-15*NH 

C2 a 1 •ESE-IS^NH^EXPC-.TS^ST ) 

C3 » 3»67E-16*NH*EXP (-•88B8*ST ) 

C4 * 1 •55E-16*NH-»EXP (-•9374*ST ) 

D1 = 1 •0E-6*(NE**»286) 

IF(D1#GT.«38 ) Dl=s*38 

A1 s 3* 16E-20*NH*T*EXP{-( 1 • 0-D 1 / 1 3 • 595 ) *ST ) ♦ ( EXP ( ( .04-01 / I 3*595 )* 
1ST)-1 .0) 

CD = 1 .35E-35»NHP*Ne/SQRT(T ) 

A( I ) s F3( 1 )*C1 
A(2) = F3(2)*C1 
CC * CD+C2+C3+C4 
W = ( 1 .0-.0736*01 >**2 
A(3) = F3(3)*CC + .014*C1# ( 1 .0-W)/W 
00 20 1=4*26 

20 Ail) - F3( I )*CC 
CC * CC-C2+A1 
W = ( 1 .0-.294-»Dl )»*2 

A(29) = F3<29)*CC + .23»C2* ( 1 • 0-W ) /W 
DO 30 1=30.44 

30 AU ) * F3< I )*CC 
CC = CC-C3 

W » ( 1 .0- *662*01 1**2 

A(45) = F3(45)*CC + 1 • 1 53*C3* ( 1 • 0-W )/W 
00 40 1=46.48 

40 A( I ) = F3( I )*CC 
CC * CC-C4 
X = *85-01 
DO 45 1=49.53 

IF(F1 ( I )-X ) 46.47.48 

46 W = (X/Fl ( I-l ) )**2 

W1 = W*F1 ( I-l )*F1 ( I-l )*(F1 ( I-l )-Fl ( I ) ) 

All) « F3(1)*CC + C4*W/\in 
K1 = 1+1 
GO TO 49 

47 A(I) = F3 ( I )* (CC+C4 ) 

K1 = 1+1 

GO TO 49 



48 A( I ) = F3( I )*(CC+C4) 

45 CONTINUE 

K1 = 54 

49 CONTINUE 

DO 50 I=K1 «58 

50 A ( I ) = F3 ( I )*CC 
C 

C *** BOUND-BOUND TRANSITIONS FOR HYDROGEN *** 

C 

51 = 1.098E“16*NH 

S( 1 ) = .4160»S1 
S(2) = .OTRl^Sl 
5(3) = •0290*51 
S(4) = .0139*51 

52 = SI *EXP (-•75*ST ) 

S(5) = 2.560*32 
S(6) = 0.447*S2 
S(7) = 0.179*S2 

S2 = S1*EXP(“.8886*ST ) 

S(8) = 7.578*52 

S(9) = 1.355*52 

52 = 51 *EXP(-.9374*ST ) 

5(10) = 16«61*52 

G(l) = 4.53E-15* (NE**.6667 ) 

G ( 1 ) = .42*G ( 1 ) 

G(2) = 4.30*G( 1 ) 

G (3 ) = 8.1 2*G ( 1 ) 

G(4 ) = 13.4*G( 1 ) 

G(5) = 1 .80*G( 1 ) 

G (6 ) » 6.44*G ( I ) 

G(7) = 11 .5*G( 1 ) 

G(8) = 2.58*G( 1 ) 

G(9) = 8.63*G( 1 ) 

G( 10 )= 3.33*G( 1 ) 

DO 300 I=3«25 

PS = 0. 

DO 200 J=l *4 

PS = PS + (.3183*S(J)*{ATAN((F1(I-1)-V0(J))/G(J))-ATAN((F1(I)- 
IVO(J) )/G(J) ) )/(Fl (I-l )-Fl ( I ) ) ) 

200 CONTINUE 

A ( I ) = A( I > + PS 
300 CONTINUE 

DO 301 1=29.44 

PS = 0. 
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DO 201 J=5»7 

PS = PS + ( J)*(ATAN( (FI ( I-l )-V0( J) )/G( J) )-ATAN( (FI ( 1 ) 
1V0(J> )/G(J) ) )/(Fl ( I-l )-Fl (I ) ) ) 

201 CONTINUE 

A( I ) = A( I ) + PS 

301 CONTINUE 

DO 302 I=45#48 

PS = 0« 

DO 202 Js5«9 

PS * PS + { .3183*S(J )* (ATAN( (FI ( I-l )-VU( J ) )/G( J) )“AIAN( (FI ( I ) 
IVO(J) >/G(J) ) )/(Fl ( 1-1 )-Fl ( I ) ) > 

202 CONTINUE 

A( I ) = Ad) + PS 

302 CONTINUE 

DO 303 1=49*58 

PS = 0* 

DO 203 J=5*10 

PS = PS + ( .SlSS^Sl J)*(ATAN( (FI ( I-l )-V0( J) >/G( J) )-ATAN( (FI d ) 
1 V0( J ) )/G< J) ) )/<Fl ( I-l )-Fl ( I ) ) ) 

203 CONTINUE 

Ad ) s Ad ) + PS 

303 CONTINUE 
990 CONTINUE 


**♦ CORRECTIONS FOR INDUCED EMISSION ♦** 

M = 58 

DO 500 I =1 »M 

Ad) = A( I )*( 1 .0-EXP(-F( I )#ST) ) 

500 CONTINUE 

return 

END 

SUBROUTINE BLBODY 

♦** CALCULATION OF BLACK0OOY FUNCTION *** 

COMMON/T 1 /T • WHO • XH2 • XHE * P • NH2 * NH « NHP • NHN * NE * NHE * NHEP » NHE2P 
COMMON /Rl/ A(58) *B(58)*F1 (58)*F3(58) *F(58) * VOdO)«Wl(58) 
REAL NH2 * NH « NHP *NHN * NE * NHE * NHEP ♦ NHE2P 
BL(X) = •15399*(6** (X+1 • ) + X*X*(X + 3. ) )*EXP(-X) 

ST = 1 57780. 0/T 
M = 58 
N = M-1 
M N1 = M+1 

DO 10 1 = 1 *N 
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Y = Fid )*ST/13*595 
IF(Y«LT*2.0) GO TO 5 

B<I> = BL(Y) +,0625*BL (2«*Y) + . 0 I 235*BL ( 3 • ♦ Y ) 

GO TO 10 
5 CONTINUE 

8(1) = 1.0 -.05133*Y*Y*Y*( 1 .0“.375*Y+.05*Y*Y ) 
10 CONTINUE 

B(M) = 1.0000 
DO 20 Isi,N 
J = Nl-I 
K = J«1 

B(J) = B(J)-B(K) 

20 CONTINUE 
RETURN 
END 


SUBROUTINE RADAT 

COMMON /CONST/ GAMA ♦ DY ( 5 1 ) ♦ Y < 5 I 
COMMON /C0NST6/ YH.YKfTB 
COMMON /DEPEND/ T < 20 ♦ 5 1 ) * P ( 20 * 5 
ll).XHH(51)* XH(5l ) *XHP(51 ) *CP<5 
COMMON /SHOCK/ OS (20)* TS(20)* 
120 ).ARFA(20 ) .CPS (20 ) 

COMMON /INF/ VINF* D INF* CP I 
COMMON /INFI/ TINF* PREI* VELOI 
COMMON /RAD/ COFl (58*51 ) * OFXP ( 
COMMON /Rl/ 0(58) *B(58 ) *F1 (58) « 
COMMON /SDFDIS/ NS ( 20 ) * NS 1 « NS2 * 
D I MENS ION V 1 ( 58 ) ♦ V2 ( 58 ) * 

1 EI2(51)*EI1(51)*E2(51)»QP1(51 

2(51 ) *EM2 (51 ) *QM1 (51 ) *W1 (58 ) *W2 ( 

real ns 

external FX *FX1 

1QW=0 


) *CC3 (20.51 ) .CCCl 

1 )*D(20*51 )*U(20*51 ).V(20*51 ) 
1 ) *XHEP(51 ) *XHE2P(51 ) *XEMIN(5 
VS(20)* US(20)* PS (20 ) ,MUS(20 


« DENI* M* ALT* 1111* CP I I 
51 ) *QFXM(51 ) »QFX(51 ) * DQY (51 ) 
F3 ( 58 ) *F (58 ) * VO ( 1 0 ) * W 1 ( 58 ) 
PSl * AS (20 ) 

C0F3 (58.51 ) * C0F2 (58*51 ) * 

) * A ( 51 ) *E3 (51 ) *QPT2 (51 ) *EJl (5 
58) *DD(50 ) *W3( 58 ) «W4 (58 ) 


K1 1 =50 

CCCl =DINF-*VINF**3 
A0=0. 26777343 
A1 =8.6347608925 
A2= 18. 05901 69730 
A3=8. 5733287401 
80=3.958469228 
B1 =21 .09965308 
82*25.63295614 


*XHE(5 
1 ) 

) * TKS ( 


1 ) *EJ2 
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63=9.5733223454 

DO 21 N=1 .50 

DD(N )=0Y(N)*23*NS (M ) 

21 CONTINUE 
YE=0.8 

YH=6.625E“27 
YK=1 .380E-16 

XX=2.*YH* (YK/YH)**4/ (3.E 10**2 ) 

M5 = 58 
EPSaO.Ol I 
DO 66 N=l .51 

301 TSS=T(M.N)*TS(M )*VINF**2/CPI 
DO 2 1=1.57 

W2 ( I )=W1 ( I+l ) 

2 CONTINUE 
W2(58)=0.01 
DO 19 1=1.58 

VI ( I )=YH*W1 ( I )/(YK*TSS)*3.E18/1 .24021 1 5E4 
V2( I )=YH*W2( 1 >/(YK*TSS)*3.E18/1 .24021 1 5E4 
W3( I )=Wl < I )*3.E18/1 .24021 1 5E4 
W4< I )=W2( I ) *3. El 8/1 .24021 1 5E4 
19 CONTINUE 

DO 67 K=1 .M5 
XXI =XX*TSS**4 

93 CALL ROMBS ( V2 (K ) . V I (K ) . FX . EPS . SUM .IERR) 

C0F2 (K.N)=XX1*SUM 
IF (IQW.EQ.O) GO TO 67 

CALL ROMBS( W4 (K ) . W3 (K ) .FXl .EPS.SUN .IERR) 

XX2 = YE*YH/(3.E1 0**2 ) 

C0F3 (K . N ) =XX2*SUN 
67 CONTINUE 
66 CONTINUE 

DO 1 N=2.51 

IF (N.GT.40) GO TO 302 

ND=N/2 

ND=ND*2 

IF (ND.NE.N) GO TO 1 

302 N1=N-1 
N2=N+l 
QPPP=0 
QMMM=0 

DO 10 K=1 .M5 
DO 9 I = 1 .N1 
J=I + 1 
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Ell ( J) = DD( I )*<C0F1 (K» I )+COFl (K* I + l ) )/2# 

9 CONTINUE 
DO 3 I = I »N 
EI2( I )=0* 

J= I 

4 IF (J.EQ.N ) GO TO 3 
EI2( 1 )=EI I ( J+1 )+EI2 ( I ) 

J = J+l 

GO TO 4 
3 CONTINUE 

EI2(N) =0«04*DD(N1 )*COFl (K*N) 

DO 5 I = I . N 

IF (EI2( I )*GT.l • > GO TO 70 

E2 ( I ) = 1 •+ (0,5772-1 *+ALOG(EI2 < I ) ) ) *E I 2 ( I ) -E 1 2 ( I >**2/2.+E12 ( I )**3/l 2 

1 • 

GO TO 5 

70 X=E12 ( I ) 

IF (X.GT.200) GO TO 152 

El 1 1 =EXP (-X )* ( AO + Al *X + A2*X**2+A3*X**3+X**4 ) / ( X* ( BO+B 1 *X-»-B2*X**2+B3 
1*X**3+X*»4 ) > 

E2 ( I )=EXP(-X )-X*El 1 1 
GO TO 5 
152 E2 ( I )=0, 

5 CONTINUE 
QPl (1 )=0, 

DO 7 I=2*N 

QPl { I )=DD( I-l )/2» (COFl (K, I-l )*C0F2(K . 1-1 )*E2 (1-1 >+COFl (K, I )*COF2 (K 
1 , I )*E2 ( 1 ) )+QPl ( I-l ) 

QPP-QPl ( I ) 

7 CONTINUE 

QPPP=QPP+QPPP 
10 CONTINUE 
QPT2 (N )=0, 

IF ( IQVK.EQ.O ) GO TO 210 
DO 100 K=1 ,M5 
A ( 1 )=0* 

DO 101 I=2iN 

A( I )=DD( I-l )/2,*(COFl (K* I-l )+COFl (K , I ) ) + A ( I - 1 ) 

AI I=A( I ) 

101 CONTINUE 

IF (AII.GT,!#) GO TO 71 

E3(N)=0.5-AI I +0,5* (-0.5772 + 1 .5-ALOG(AI I > )*AI 1**2 + A I l**3/6» 

GO TO 72 

71 X-AII 


IF (X.GT*200) GO TO 151 

Ell 1 =EXP < -X ) ♦ ( AO+ A I *X+A2*X**2+A3*X**3+X**4 ) / ( X* < BO+Bl #X+B2*X**2+83 
1*X**3+X**4 ) ) 

E222 = EXP(“X)~X*E1 1 1 
E3(N)=(EXP(-X)-X*E222)/2« 

GO TO 72 
151 E3(N)=0* 

72 CONTINUE 

QPT2(N)=E3(N)*COF3(K*N)+QPT2 (N) 

100 CONTINUE 
210 QMMM=0* 

IF (N.EQ.51 ) GO TO 190 
DO 112 K=1 f MS 
DO 113 I=N1 *50 
J=I + 1 

EJl ( J > = DDCI )/2**(C0Fl (K* I )+COFl (K* I + l ) > 

1 13 CONTINUE 

EJl (N1 )=0,040»DD(N1 >*C0F1 (K,N) 

EJ2(N )=0* 

DO 114 I=N2*51 
IF {N2*EQ.52> GO TO 114 
EJ2( I )=EJ1 ( 1 )+EJ2 (1-1 ) 

114 CONTINUE 
EJ2(N)=EJ1 (N1 } 

DO 115 I=N *51 

IF (N2.EQ.52) GO TO 115 
IF (EJ2 ( I ) .GT* 1 • ) GO TO 74 

EM2( I )=1 •+< .5772-1 .+AL0G(EJ2( I ) ) )*EJ2( I )-EJ2 ( I )**2/2.+EJ2( I )**3/12 

1 • 

GO TO 115 
74 X=EJ2 ( I ) 

IF (X.GT.200) GO TO 150 

El 1 1 =EXP(-X )* (A0+A1*X+A2*X**2+A3*X**3+X**4 ) / ( X* ( BO+B 1 *X+B2*X**2+83 
1*X*#3+X**4 ) ) 

EM2( I )=EXP(-X)-X*E1 1 1 
GO TO 115 
150 EM2( I )=0. 

GO TO 115 

1 15 CONTINUE 
QMl (N-1 )s=0. 

D0116 I=N2*51 

IF {N2.EQ.52) GO TO 116 
IP=I-1 

QMl ( IP)=DD( I-l )/2*(C0Fl (K* I-l )#C0F2<K* I -I )*EM2( I-l )+COFl (K* I )*C0F2 


I-* 

to 

-j 
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1 <K * I )*EM2 ( I ) ) + QMl ( IP-1 ) 

160 QMMsQMl ( IP ) 

116 CONTINUE 
1 03 QMMM=QMM+QMMM 
112 CONTINUE 

190 QFXP(N)=(QPPP+QPT2(N) )*6.2831/CCC1 
00 QFXM (N ) =QMMM*6«2831 8/CCC 1 

QFX (N ) =QFXP (N >-QFXM (N ) 

1 CONTINUE 

DO 303 N-3»39,2 

QFXP(N )- (QFXP(N+1 )+QFXP(N-l ) )/2 
QFXM (N ) = (QFXM (N+1 )+QFXM (N-1 ) )/2 
QFX(N >= (QFX (N+1 )+QFX(N-l ) )/2 
303 CONTINUE 

QFXM ( 1 ) sQFXM (2 ) 

QFXP( I )=0.8*5.668E-5* (TB #*4 )/CCCl 
QFX ( 1 ) =QFXP ( 1 ) -QFXM ( 1 ) 

DQY( 1 )= (QFX(2)-QFX( 1 ) )/OY( 1 )/NS(M) 

DQY(51 )=“(QFX(50)-OFX(51 ) ) /DY ( 50 ) /NS ( M ) 

DO 888 I=2»K1 1 

AAl=DY( I-l )/(DY ( I )*(DY ( I-l )+DY( I ) ) ) 

AA2=DY (I )/ ( DY ( 1 -I )* ( DY ( I -1 )+DY ( I ) ) ) 

AA3=(0Y( I )-DY( I-l ) )/(DY( I ) *DY ( I-l ) ) 

DQY( I )=AA1*QFX ( I+l )-AA2*QFX( I-l >+AA3*QFX( I ) 
889 DQY( I )=DQY( I )/NS(M) 

888 CONTINUE 
RETURN 
END 

FUNCTION FX(X) 

FX=X**3/(EXP(X)-1 • ) 

RETURN 

END 

FUNCTION FXl (X ) 

COMMON /C0NST6/ YHfYK*TB 
FX1=X»*3/(EXP(YH*X/(YK*TB> )-l ) 

RETURN 

END 


SUBROUTINE ENERGY 

COMMON / BODY/ 2TA(20)« CK(20)* X(20)» DX • R ( 20 ) » BATA ( 20 > • MUREF. EPS 
COMMON /CONST/ GAMA . DY ( 5 1 ) » Y ( 5 1 ) t CC3 ( 20 « 5 1 ) »CCC1«IC0NT 
COMMON /CONSTl/ AL * AK . CT * K • K 1 1 « MBP * MA 1 * M A2 

COMMON /C0NST3/ AH2 ( 3 ♦ 7 ) » AHP ( 3 , 7 ) » AH2P ( 3 ♦ 7 ) t AHE ( 3 , 7 ) ♦ AHEP ( 3 * 7 ) ♦ 
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1 AE(3*7) tAH(3»7)fKE 

COMMON /C0NST6/ YH*YK*TB 

COMMON /DEF/ DVY(51)» DK(51)« DMU(51)« DPY(51)* DUY(51>* DPI (51) 
COMMON /DEPEND/ T(20,5n»P(20»51)tD(20,51 )»U(20*51 )iV( 20*51 ) *XHE(5 
1 1 ) «XHH(51 ) * XH(51 ) *XHP(51 ) *CP(51 ) .XHEP(5I ) *XHE2P(51 ) *XEMIN(51 ) 

1 »THP1 (51 ) *TH (51 ) 

COMMON /ENTH/ THH2* THH* THHE* THHP* THE 

common /heat/ qcon*difu* rat ,QT0TL 

COMMON /INF/ VINF* D INF* CP I 

COMMON /INFI/ TINF* PREI* VELOI* DENI* M* ALT* Mil* CPII*XH2I 
1 » NORAD *NOPRE 

COMMON /PERC/ V I ( 20 ) ,P I NF ( 20 ) ♦ ETHP(20), D I ( 20 ) ♦ CPL ( 20 ) . 001 (20)* 

1 PINFN(20) 

COMMON /RAD/ COFl (58*51)* QFXP ( 51) ♦ QFXM ( 5 1 ) ♦ QFX ( 5 1 ) * DQY ( 5 1 ) 

COMMON /SOFDIS/ NS ( 20 ) * NS 1 * NS2 * PS 1 ♦ AS ( 20 ) 

COMMON /SHOCK/ DSN ( 20 ) * TSN ( 20 ) ♦ VSN ( 20 ) ♦ USN (20 ) * PSN ( 20 ) * MUSN(20)* 

1 TKSN (20 ) . ARFA (20 ) ,CPSN (20 ) *THS (20 ) 

COMMON /SPHT/ CH2 ( 5 1 ) * CH ( 5 1 ) * CHP ( 5 1 ) * CHE ( 5 1 ) 

COMMON /TRANS/ TK ( 5 1 ) ♦ MU ( 5 1 ) • TKD I 

DIMENSION A2 (51 ) » A3 (51 ) * A4 (51 ) * AN(5l ) *BN( 51 ) *CN (51 ) *DN (51 ) * 

1E(51)*F(51) *DPX(51 ) *DNS(20) *DTS(20) *DPS(51 )*A1 (51 ) 

2 *DHP(51)* DHE(5l)i DH(51)* DH2(51)* PR ( 5 1 ) * DPR ( 5 1 ) *DE(51) 

3 *CH2P(51 ) *CE(51 ) * DPCI (51 )*DTH(51 ) • PC I ( 5 1 ) * DTHS ( 20 ) 

REAL MUHH* MUHE * MUH t MUD . MSD ♦ MU* MAC I » MUREF » MUSN • NS ♦ NS 1 * NS2 

IF (N0RAD*NE*1) go TO 31 

DO 32 N=1 *51 

QFXP(N ) =0 

QFXM (N ) =0 

0FX(N)=0# 

DQY(N)=0 
32 CONTINUE 
31 HX2I=XH2I 
V(M*51 ) = 1 • 

U(M*51 )-l • 

D(M*51 ) = 1 • 

P(M*51 ) = 1 • 

TH(5l ) = 1 • 

T(M*51 ) = 1 • 

T(M*1 )*TB*CPI/ (VINF*VINF*TSN(M) ) 

PRS=CPSN ( M ) *MUSN ( M ) /TKSN ( M ) 

DO 131 N=l* 51 

AAA=XHH (N)*2+XH(N)+XHP(N)+XHE(N)*4+XEMIN(N)*0.5A86E-3 

CH2(N)=XHH(N)»2/AAA 

CH(N)= XH(N)/AAA 
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CHP(N ) =XHP(N>/AAA 
CHE(N ) =XHE(N>/AAA*4. 
CE(N)=XEMIN(N)*0,5486E“3/AAA 
PR(N)=CP(N)*MU(N)/TK (N) 

CH2P(N )=0« 

131 continue 

DO 2 N1 =1 *49 
N=51-N1 

AA4 = DY(N)*( DY (N-1 )+DY (N ) ) 

AAl-OY (N-1 )/AA4 

AA2 = -DY (N)/(DY(N-1 )* (DY(N-*1 )+DY(N) ) ) 

AA3 = +(DY(N)-DY(N-'l ) ) / ( DY ( N ) *DY ( N- 1 ) ) 
DMU<N)=AAl*MU(N+l ) + AA2*MU ( N- 1 ) A A3*MU ( N ) 
DK(N)=AA1*TK(N+1 )+AA2»TK(N-1 )+AA3*TK(N) 
DPY(N)=AA1*P(M,N+1 )+AA2*P(M,N-1 )+AA3*P(M*N) 
DUY(N)=AA1*U(M,N+1 )+AA2*U(M«N-1 )+AA3*U(M»N) 
DVY(N)=AA1*V(M,N+1 >+AA2*V (M,N-1 }+AA3*V(M*N) 
DHP(N)=AA1*CHP(N+1 )+AA2* CHP(N“1 )+AA3*CHP(N> 
DHE (N > =AA1 *CHE (N+1 )+AA2* CHE(N“1 )+AA3*CHE(N) 
DE (N)=AA1*CE (N+1 )+AA2* CE (N-1)+AA3*CE (N> 
DHa(N)=AAl*CH2(N+l )+AA2* CH2<N“1 )+AA3#CH2(N) 
DH (N)=AA1*CH (N+1 )+AA2* CH (N“1)+AA3*CH <N) 
DPR(N)=AA1*PR (N+1)+AA2* PR (N-1 )+AA3*PR (N) 
2 CONTINUE 

DHE { 1 )= ( (CHE (2 ) -CHE ( 1 ) )/DY ( 1 ) +DHE(2 ) )/2 
DHP( 1 )= ( (CHP(2)-CHP( 1 ) )/DY( 1 ) +DHP ( 2 ) )/2 
DH2 (!) = (( CH2 ( 2 ) -CH2 ( 1 ) ) /OY ( 1 ) +OH2 ( 2 ) ) /2 
OH (1>=((CH (2)-CH (1))/DY(1> +DH(2))/2 
DE (1)=((CE (2)-CE (1))/DY(1) +DE(2))/2 
OPR( I >= <PR(2)-PR< 1 ) )/DY( 1 ) 

DPR (51 )=0. 

DMU( 1 )= (MU(2 )-MU( 1 > >/DY( 1 ) 

DK ( 1 )= ( TK (2 l-TK ( 1 ) )/DY ( 1 ) 

DP Y (1 )=(P(M*2)-P(Mtl )) /D Y ( 1 ) 

DUY (1 )~(U(M*2)-U(Mfl )) /DY ( 1 ) 

DK(51 )=-(TK(50)-TK(5I ) >/DY(50) 

DPY(51 ) =-(P(M,50 )-P(M»51 ) )/DY(50 ) 

DUY (51 )=- (UCMfSO )-U(M«51 ) ) /DY ( 50 ) 

DVY( 1 >=(V(M,2)-V(M,1 ) )/DY( 1 ) 

DMU(51 )s (MU<5! )-MU(50 ) )/DY(50 ) 

DVY(51 )=(V(M»51 )-V(M,50) )/DY(50) 

KE=1 

DO 101 N=UKll 

TEMP=T(M.N)*TSN(M )*V I NF**2/CP I 
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call ENTHALP (TEMPtCH2(N) *CH(N) » CHE ( N > f CH2P ( N ) *CE(N) *CHP(N) tCPC* 

1 ENTHC) 

KE = 2 

IF (N.NE. 1 > GO TO 106 

THM )*(ENTHC/<VINF»*2 ) )/THS(M) 

106 AND0= VINF*»2 

THH2=THH2/ ( 2*AND0 ) 

THHE=THHE/ ( 4*AND0 ) 

THE=THE/(0*5486E“3*AND0 ) 

THH=THH/AND0 

THHP=THHP/ANDO 

SGM = DH2 ( N ) ♦THM2+DH < N ) *THH+DHE ( N ) *THHE+DHP ( N ) *THHP +OE ( N ) #THE 
SGM=SGM/ (2*390 lE-e ) 

IF (N *NE* 1 ) GO TO 54 

DIF=-1 • 1*MU<N)*MUSN(M)»SGM/NSCM ) /PR ( N ) 

54 PCI=MU(N)/PR(N)*0*1*SGM 

PC2=MUSN ( M ) *USN ( M ) **2*CK ( M ) #MU ( N ) *U ( M ♦ N ) *#2/ ( 1 +NS ( M ) »Y ( N ) *CK ( M ) ) 
PC3sMU ( N ) *U { M ♦ N ) *USN ( M ) ** 2 * ( PRS*PR ( N ) -1 > *DUY ( N ) /PR ( N ) 

PC I ( N ) = MUSN ( M ) / ( NS ( M ) *PRS ) * ( PC 1 +PC3 ) -PC2 

101 CONTINUE 
KK=K1 1“1 

DO 102 N=2*KK 

AA4=DY ( N ) * ( DY ( N- 1 ) +DY ( N ) ) 

AA1=DY(N-1 )/AA4 

AA2=-0Y (N )/<DY(N-l )*(DY(N-1 )+DY(N) ) ) 

AA3=+(DY(N)~DY (N-1 ) )/(DY(N)*DY(N-l ) ) 

DPCI (N)=AA1 *PCI (N+1 )+AA2*PCI (N-1 )+PCI (N)*AA3 

102 CONTINUE 

DPCI (K1 1 )sDPCI (KK ) 

DPC I ( 1 ) = ( PC I ( 2 > -PC 1(1)) /DY ( 1 ) 

IF (M.GT* 1 ) GO TO 5 
NS1=NS( 1 ) 

PS1=PINFN(M) +( 1 ,-l ./DSNIM) ) 

DO 6 N=1 *K1 1 

A81-DMU (N )/MU(N )-DPR (N ) /PR ( N ) +2*NSl / ( 1+NSl #Y (N ) ) 

A 1 ( N > = AB 1 -NS 1 *DSn < M ) *PRS* VSN ( M ) *D ( M » N ) *PR ( N ) * V ( M » N ) / (EPS»*2»MUSN ( M 
1 )*MU(N) ) 

A4 (N ) = 0 
A2 (N )=0 • 

AB2=PRS»NS 1 **2*PR ( N ) / ( MUSN ( M ) *THS ( M ) *MU ( N ) ) 

A3(N)sA82*(DPCI (N )/NSl+2*PCI ( N ) / ( 1 +NS 1 *Y ( N ) ) ) 

A3(N)=A3(N>-(A82/EPS»*2)*(DQY(N)+QFX(N)*(2/(1+ NSl ♦Y(N))) ) 

A3(N)=sA3(N)+DPY(N)*VS N( 1 )*V(1 *N)*PSN( 1 ) * AB2/ ( EPS**2*NS 1 ) 

6 CONTINUE 
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GO TO 10 
5 DO 11 N=1 *K1 1 

DPX(N)= (P(M,N)-P(M-i ,N J )/0X 

11 CONTINUE 

IF (M,EO*MBP) GO TO 59 

DTHS(M )= (THS(M+1 )-THS<M-l ) )/<2*0X) 

IF (1111 .EQ#! ) GO TO 50 
DNS(M)= (NS(M+1 )-NS(M-l ) )/(2**DX) 

GO TO 52 

59 IF (Illl.EQ.l) GO TO 50 

DNS(MBP)=(AS(MBP)-AS(MA2) l/DX 
DTHS(MBP) = (THS(MBP)“THS(MA2 ) )/0X 
GO TO 52 
50 DNS(M)=0 

IF (M#NE.MBP > GO TO 52 
DTHS ( MBP ) = ( THS ( MBP ) -THS ( MA2 ) )/DX 
52 DO 12 N*1 fKl 1 

TCI=1+ Y{N)»CK(M)#NS(M) 

TC2*EPS**2*MUSN ( M )*MU ( N ) 

TC3*C0S (ZTA (M) >/ ( R ( M ) +NS ( M ) * Y ( N ) *COS ( ZTA ( M ) > > 

AB3 = DMU (N)/MU(N)-DPR(N)/PR(N)+NS(M)*(CK(M >/TCl+TC3) 

AB4= (DNS<M)*USN(M )*Y(N)*U(M*N)*D (MtN)/TCl -VSN(M)*D (MtN)*V (M*N> ) 
A1 (N)=AB3+DSN(M)*PRS*PR(N)*NS(M )/TC2*AB4 
AB5=PRS»PR ( N ) *NS ( M ) **2/ ( MUSN < M ) *MU ( N ) *THS < M ) ) 

A3(N)=AB5*(DPCI (N )/NS ( M )+ (CK <M )/TCl +TC3 )*PC 1 (N) ) 

A3(N)=A3(N)- (AB5/EPS**2)* (DQY(N )+QFX(N>*(CK<N)/TCl+TC3) ) 
A3(N)=A3(N)+DPY (N)*AB5*VSN <M)*V (M,N)/NS(M )*PSN(M ) /EPS*»2 
66 A4 ( N ) =PRS*NS ( M > **2»DSN ( M ) »USN ( M ) *PR ( N ) *U ( M f N ) *D ( M • N ) / < TC 1 *TC2 > 
A4(N)5-A4(N) 

A2 <N ) s A4 ( N ) ♦DTHS ( M ) /THS ( M ) 

12 CONTINUE 

10 DO 20 N=2»K1 1 

AN(N)= (2*+Al (N)#DY(N-1 ) ) / ( DY (N ) +DY ( N- 1 ) )/OY(N) 

BN(N)=(2.-A1 (N)*(OY (N>-OVfN~l ) ) ) / ( OY ( N ) *DY ( N- 1 ) ) -A4 < N ) /DX-A2 ( N ) 
BN(N )«-BN(N ) 

CN(N)s(2#-Al (N)*DY(N) )/(DYCN-l ) * ( DY < N > +DY ( N- 1 ) ) ) 

IF (M.GT.l ) GO TO 18 
DN(N)=A3(N) 

DN (N ) =-DN(N ) 

GO TO 20 

18 DN(N)=A4<N)/DX*THP1 (N) *(-l)+A3(N) 

DN(N)=-DN(N) 

20 CONTINUE 
Ed )=0« 
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F(1 )sTH(l ) 

DO 22 N=2»K1 I 

E(N>=-AN(N)/(BN(N)+CN(N)*E(N“1 ) ) 

F(N>= (DN(N)-CN(N)*F(N“1 ) ) / ( BN ( N )+CN ( N ) #E (N-l ) ) 

22 CONTINUE 

DO 24 Nl=l »49 
N-51-NI 

TH(N)=E(N)*TH(N+1 >+F(N) 

24 CONTINUE 
YH2P=0* 

DO 108 N4=2»50 
N»52“N4 

RHO=D(M*N)*DSN(M)*DINF»1000 

145 TSS2sT(M»N)#TSn(M )*V INF*»2/CP1 

146 call NOBDEN (TSS2* RHO » H2 * H « HPLUS« HE * HE PL US ♦ HE2PLUS*EM INUS< 

1 XH2 I ) 

call VOLMAS (H2*H» HPLUS«HE» HEPLUS*HE2PLUS.EMINUS« YH2* YHt YHP « 
1 YHE* YHEPiYHE2Pi YE) 

CALL ENTHALP ( TSS2 * YH2 * YH* YHE * YH2P » YE * YHP» CPC * ENTHC ) 

THO= ( ENTHC/V I NF**2+U ( M * N ) *»2/2*USN ( M ) **2 ) 

137 DTHO=THO-TH(N)*THS(M ) 

TSS3= TSN(M)*VINF»*2/CP1*(T (M*N +l)-0*05) 

110 call NOBDEN (TSS3* RHO* H2 ♦ H * HPLUS * HE • HEPLUS* HE2PLUS*EM INUS* 

1 XH2 I ) 

CALL VOLMAS (H2*H« HPLUStHE* HEPLUS ♦ HE2PLUS * EM I NUS * YH2 * YH* YHP* 
1 YHE* YHEP*YHE2P* YE) 

CALL ENTHALP < TSS3 * YH2 * YH « YHE • YH2P « YE * YHP *CPC * ENTHC ) 

THl = < ENTHC/V I NF**2 + U ( M ♦ N ) **2/2*USN ( M ) **2 ) 

136 DTH1=TH1-TH(N)*THS(M ) 

111 IF (ABS(DTH1/Th1 ) .LT,0,01 ) GO TO 109 
SL=(DTH1-DTH0)/<TSS3-TSS2 ) 

TSS2=TSS3 

DTH0=DTH1 
TSS3=TSS2-DTH0/SL 
GO TO 110 

109 T(M*N)*TSS3»CPl/tVINF*VINF»TSNlM> ) 

108 CONTINUE 

DTl = <T(Mt2)-T(M* 1 ) )/DY(l ) 

DT2=DT1 *TSN<M )/NS (M ) 

QCON=TK < I )*TKSN (M )* DT2*EPS**2* ( -1 ) 

DIFU=DIF*EPS**2 
QT0TL=-1 * (QCON+DIFU ) 

RETURN 

END 
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SUBROUTINE VOLMAS ( H2 * H * HPL US « HE * HEPL US t HE2PL US , EM I NUS * CH2 * CH t CHP 
1 *CHE*CHEP» CHE2P,CE ) 

AA 1 =H2+H+HPLUS+HE+HEPLUS+HE2PLUS+EM I NUS 

XH2sH2/AA1 

XH=H/AA1 

XHP=HPLUS/AA1 

XHE=HE/AA1 

XHEP=HEPLUS/AAl 

XE=EMINUS/AA1 

AA2=XH2*2+XHE*4+XH+XHP+XHEP*4+XE*0.5486E-3 

CH2=XH2*2/AA2 

CH=XH/AA2 

CHP=XHP/AA2 

CHE=XHE*4/AA2 

CHEP=XHEP*4/AA2 

CE=XE*O*5406E-3/AA2 

RETURN 

END 


SUBROUTINE MOMENTM 

COMMON / BODY/ ZTA(20>* CK(20)t X(20)f DX • R ( 20 ) ♦ BATA ( 20 > * MUREF ♦ EPS 
COMMON /CONST/ GAMA * DY ( 5 1 ) ♦ Y ( 5 1 ) » CC3 ( 20 » 5 1 ) «CCC1*IC0NT 
COMMON /CONSTl/ AL < AK * CT . K t K 1 1 • MBP » M A 1 » M A2 

COMMON /DEF/ DVY(51)» DK(51)t DMU(51)* DPY(51)* DUY<51)« DPI (51) 
COMMON /DEPEND/ T(20*51 )»P(20*51 )«D(20,5l )*U(20«51 )»V(20.51 ) *XHE(5 
11)* XHH (51)* XH ( 5 1 ) * XHP (51 ) * CP ( 5 1 ) 

COMMON /INFI/ TINF* PRE I * VELOI ♦ DENI* M» ALT* Mil* CPI I 
COMMON /INF/ VINF* D INF* CP I 

COMMON /PERC/ V I ( 20 ) *P I NF ( 20 ) * ETHP(20)* D I ( 20 ) ♦ CPL ( 20 ) * Q0l(20)» 

1 PINFN(20) 

COMMON /SDFDIS/ NS ( 20 ) * NS 1 * NS2 * PS 1 * AS ( 20 ) ,NSS2 
COMMON /TRANS/ TK ( 5 1 ) * MU ( 5 1 ) 

COMMON /SHOCK/ DS(20)* TS(2Q)* VS(20)* US(20)* PS ( 20 ) * MUS ( 20 ) « TKS ( 

120 ) ,ARFA(20 ) 

DIMENSION PI (51 )*P2 (51 )*U1 (51 )*E1 (51 ) *F1 (51 )*B1 (51 )*B2(51 )*B3(51 )• 
1B4 (51 ) *DPX(51 ) * ANl (51 ) *BN1 (51 ) ,CN1 (51 ) »DN1 (51 ) * DNS (20) * 

2 DPS(20)* DVS(20) ,DVX(5l ) *DUS(20) 

REAL MUREF* MUHH * MUHE • MAC I » MU * MUS • MUH * NSS 1 , NSS2 * NSS3 
real NS*NS1 *NS2* INTGA* INTGB* intgc* INTGD* INTGE* INTGF 

IF (M.GT. 1 )GO TO I 

PI (51 ) = 1 • 

PS1=PINFN(M) +( 1 .-1 ./DS (M ) ) 

I F (Mil .GT* 1 ) GO TO 2 



NS2 = 0 
NSl =NS< 1 ) 

PS2=-( 1 •“! ./DS(M) > 

P2 (51 )=0. 

US 1 = I • 

DO 5 N1 =l *50 
N=52“N1 

P2 (N-1 )=-DY(N-l )* (DS (M ) *US 1 **2*NS 1 »D ( M *N ) *U ( I * N ) **2/ ( PS 1 ♦ ( I *+Y(N)* 
INSl ) ) )+P2 (N) 

PI (N-1 ) = 1 • 

5 CONTINUE 
GO TO 24 

2 NS2* ( AS (3 )-AS ( 1 > )/(4.*DX**2 ) 

NS1=NS( 1 ) 

USl = 1 •-2«»NS2/ ( 1 *+NSl )* ( 1 .-1 ./DS ( 1 > ) 

PS2= ( 1 •-( 1 */DS( 1 )>)*(! .-2.»NS2/ ( 1 .+NS1 ) )**2* (“1 • ) 

P2(51 )=0. 

DO 23 Nl=l »50 
N=52-N1 

Cl 1=DS(M)*US1**2*NS1*D(M.N)*U(M»N)**2/(PS1*( 1 .+Y(N)*NS1 ) ) 

Cl 2=2*DS (M)*USl *NS2*VS (M )*D (M*N)*U(M »N)*Y (N >*DVY (N )/ (PSl * ( I .+Y(N)* 
INSl ) ) 

C13=PS2*DS( 1 > *D( I .N)*VS( I >**2*V( 1 .N)/ (PS 1**2 )*DVY(N) 

P2(N-1 )=-DY(N-l )*(C1 1+C12+C13)+P2(N) 

23 CONTINUE 

DO 99 N1 =1 iSO 

N=52-N1 

I=N-1 

PI (N-1 ) = (-VS ( 1 )**2*DS ( I )*D( 1 * I )*V ( 1 » I )*DVY( I )/PSl )*(-DY( I ) )+Pl ( N) 
99 CONTINUE 

24 DO 7 N=1 *K1 1 

B1 (N)=DMU(N )/MU(N )+2.*NSl/( 1 *+NSl*Y <N) )-NSl*DS ( 1 )*VS( 1 )*D ( 1 *N)*V( 1 
1 * N ) / ( EPS**2*MUS ( 1 ) *MU ( N ) > 

7 CONTINUE 

DO 8 N=1 *Kl 1 

BBl=DMU(N)/MU(N)+NSl/( 1+NS1*Y(N) )*2. 

BB2=DS(M )*NS1 *USl *MU ( N ) *D ( M ♦ N ) / ( EPS**2*MUS ( 1 )*MU(N) ) 

BB3=NS1 *DS( 1 )*VS( 1 )*D ( 1 *N)*V ( 1 » N ) / ( EPS**2*MUS ( 1 )*MU(N ) ) 
B2(N)=-(NS1/(1+NS1*Y(N) ) > * ( BBl +BB2+BB3 ) 

8 CONTINUE 

DO 9 N= 1 fKl 1 

H DPI (N)=(P1 (N+1 )-Pl ( N })/DY(N) 

9 CONTINUE 

DO 10 N=1 *K1 1 
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CCl =P2(N)+PS2*Pl (N)/PS1~NS2*Y(N)*DP1 (N)/NS1 

CC2=-2.*PS1*NS1**2/(EPS**2*MUS( 1 )*( 1 .+NSl»Y(N) >*US1*MU(N) ) 

B3 (N )=CC2*CC1 

10 CONTINUE 

DO I I N = 2*K11 

ANl (N) = (2*+Bl (N)*DY(N-1 ) ) / ( DY ( N ) * ( DY ( N ) +D Y ( N~ I ) ) ) 

BNl (N)s-(2«-Bl (N)*(DY(N)-DY(N-l ) ) ) / ( DY ( N 1 *DY { N- 1 ) )+B2(N) 

CNl (N)= (2«-Bl (N)*DY(N) )/(DY(N-l ) * ( D Y ( N ) +D Y ( N- 1 > ) ) 

DNl (N ) =B3 <N ) 

DNl (N ) =-DNl (N ) 

11 CONTINUE 
El ( 1 )=0. 

FI ( 1 )=0# 

DO 12 N-2*K11 

El (N )=-ANl (N)/ (BNl (N )+CNl ( N ) *E 1 (N-1 ) ) 

FI <N)= (DNl (N)“CN1 (N)*F1 (N-1 ) )/(BNl (N)+CN1 (N)*E1 (N-1 ) ) 

12 CONTINUE 
U(M,5l ) =1 • 

U (Mf 1 )=0« 

DO 13 Nl=l »49 
N=51-N1 

U(M,N) = -AN1 (N)/(8N1 (N)+CNl (N)*E1 (N-1 ) )*U(M*N+1 )+(DNl (N)-CNl (N)*F1 ( 
lN-1 ) )/(BNl (N)+CN1 (N)*E1 (N-1 ) ) 

13 CONTINUE 
I NTGA=0* 

INTGB=0# 

DO 14 N=1 «K1 I 

INTGA=DY(N)/2.*(D(M*N+1 )*U(MtN+l ) +D ( M *N ) *U < M » N ) )+INTGA 
INTGB=DY(N)/2**(D(M«N+1 )*U(MtN+l )*Y(N+1 ) +D ( M * N ) *U ( M » N )*Y ( N ) )+INTGB 

14 CONTINUE 
VSl 1 =VS (M ) 

All=(VSll +2.*US1*INTGB ) 

B11=(2.*VS11 +2** INTGA*US1 ) 

NS1= (-B1 1+(B1 1**2-4.*A1 1*VS1 1 ) **0 . 5 > / ( 2 • ♦ A 1 1 ) 

NS( 1 )=NS1+NS2*BATA( 1 )**2 
DA = 1 • 

DO 330 N=2*50 
DB= ( 1 •+NS1*Y (N) ) 

DCs=(OB*-D( 1 *N)*U( 1 «N)+DA*D( 1 iN-1 )*U( 1 *N-1 ) )/2**DY(N-l ) 

V( 1 ♦Nisi ./(DB**2»VS(M)*D(MfN) )♦ ( DA»*2*VS ( M ) *D ( M ♦ N- 1 )*V( 1 *N-1 )-2. 

1*NS1*US1*DC ) 

DA = DB 

330 CONTINUE 

DO 26 N=1 fK 
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P( 1 .N)=P1 (N )+P2 (N)*BATA( 1 )»*2 
26 CONTINUE 
GO TO 28 

1 IF Mill *GT.l ) GO TO 40 
DNS(M)=0 
GO TO 42 

40 IF (M.EQ.M8P) GO TO 33 

DNS(M)= (AS(M+1 )-AS(M-l ) )/(2**DX) 

GO TO 42 

33 DNS(MBP)=: (AS(MBP)-AS (MA2 ) )/DX 
42 DO 44 N=1 »K1 1 

01 1 =DMU ( N ) /MU { N > +NS < M ) *CK ( M ) / ( 1 • +NS < M ) *Y ( N ) *CK { M ) ) +NS ( M ) *COS ( ZT A ( M 
1 ) >/(R(M>+NS(M)*Y{N)*COS(ZTA(M) ) ) 

0 1 2=NS ( M ) *DS ( M > *US ( M ) *DNS ( M ) *0 < M » N ) *U ( M t N > * Y ( N ) / ( EPS**2*MUS ( M ) * { 1 . 
1+NS(M)*Y(N)*CK{M) )*MU(N) ) 

B1 (N)=D1 1+012-NS(M)»DS(M)»VS(M)*0{M*N)*V(M*N)/ (EPS*»2*MUS (M )*MU (N ) 
1 ) 

44 CONTINUE 

DUS(M)= {US(M)-US(M-l ) >/DX 
DO 45 N=1 ,<l I 

D2 1 =-CK <M )*NS (M >*DMU (N )/ ( ( 1 *+NS < M ) *Y(N)*CK (M ) ) *MU < N ) )-CK ( M)**2*NS ( 
lM)*-»2/< 1 •+NS<M>»Y(N>*CK(M) )**2 

D22 = -C0S (ZTA CM ) )#NS ( M ) **2*CK (M )/ (R(M )+NS ( M ) *Y ( N )*COS ( ZTA ( M ) ))/(!• + 
INSCM )»CK(M)*Y (N) ) 

D23 = -*DS C M ) *NS ( M ) **2*DUS ( M ) *U ( M « N ) *0 ( M ♦ N ) / ( EPS*#2*MUS ( M ) * ( 1 +NS CM ) *Y 
1 (N )*CK(M) )*MU(N ) ) 

D24=-NS (M)*#2*0S(M)*VS(M )*CK (M)*D(M,N)*V(M*N)/ (EPS**2*MUS ( M )* ( 1 .+ 
1NS(M)*Y(N)*CKCM) )»MU(N) ) 

B2 (N )=D21 +D22+D23+D24 

45 CONTINUE 

DO 46 N=l *Kl 1 

DPX(N)= (P(M*N)-P(M-1 ,N) )/DX 

46 continue 

DPS(M)= (PS(M)-PS(M-l ) )/0X 
DNS ( M ) = ( NS C M ) -NS ( M- 1 ) ) /DX 
DPYCK ) = -(P(M«K-l )-P(M*K) )/DY (K-1 ) 

DO 48 N=1 *K1 1 

D3 1 = { DPX ( N ) +DPS ( M ) *P CM ♦ N ) /PS ( M ) -DNS ( M ) #Y ( N ) *DP Y ( N ) /NS CM)) 

D32=~PS (M )*NS CM )*#2/ (EPS»»2*MUS CM )* ( 1 •+NS (M )*Y (N )*CK (M ) )»MUCN >*US ( 
IM) ) 

B3(N)=sD31*D32 
48 CONTINUE 

DO 49 N=1 tKl 1 

B4 (N )=DS(M)*US (M)#NS (M )**2*D (M.N )#U(M*N) / (EPS**2*MUS (M )* ( 1 • 


1 



I+NS(M)#Y(N)*CK(M) )*MU(N) ) 

B4(N)=-B4(N) 

49 CONTINUE 

DO 50 N=2*K1 I 

ANl (N)= (2*+Bl (N)*DY(N-1 ) ) / ( DY ( N ) * ( 0 Y ( N ) +D Y ( N- 1 ) ) ) 
w BNl (N ) =- (2* -B1 (N)*(DY(N)-DY ( N- 1 ) ) )/ ( DY ( N ) *DY ( N-1 ) ) +B2 (N)+B4 (N)/DX 

“ CNl (N) = (N)*DY (N) )/(DY(N-l ) * ( D Y ( N ) +D Y ( N~ 1 ) ) ) 

DNl (N)=B3(N)-B4 (N)/DX*U(M-1 tN) 

DNl (N ) =-DNl (N ) 

50 CONTINUE 
El (1 )=0. 

FI (1 )=sO. 

DO 52 N=2.K11 

El (N)=-AN1 (N)/(8N1 (N)+CN1 (N)*E1 (N-1 ) ) 

FI (N)a (DNl (N)-CNl ( N ) *F 1 ( N- 1 ) )/<3Nl (N)+CN1 (N)*E1 (N-1 ) ) 

52 CONTINUE 
U(M*Kl = I 
DO 54 Nl = 1 » 49 
N-51-N1 

U(M,N)=-AN1 (N)/(BN1 (N>+CN1 (N)*E1 (N-1 ) )*U(M,N+1 )+(DNl (N)-CNl (N)*F1 ( 
lN-1 > )/(BNl (N)+CN1 (N)*E1 (N-1 ) ) 

54 CONTINUE 
I NTGDsO • 

INTGE=0* 

ANTGM^O* 

ANTGN-0 • 

DO 56 N=1 *K1 1 

INTGD = DY (N)/2.*(D(M*N+1 )#U(M*N+1 ) +D ( M » N ) *U ( M ♦ N ) )+lNTGD 
I NTGEsDY (N)/2 t* (D (M «N+1 )*U (M t N+1 )*Y (N+1 ) +D ( M » N ) *U ( M * N ) *Y ( N ) ) + INTGE 
ANTGM = DY(N)*(0(M-1 iN+1 )*U(M-1 *N+1 )*Y(N+1 )+D(M-1 «N)^*-U(M-1 «N)*Y(N ) )/ 
1 2 • +ANT G M 

ANTGN=DY (N)»(D(M-1 ,N+1 )»U(M-1 »N+1 >+D(M-1 *N)»U(M-1 »N) )/2*+ANTGN 
56 CONTINUE 

CCl = INTGE*C0S (ZTA (Ml )*DS ( M )*US (M ) 

CC2= INTGD^R (M1*DS (M )*US(M ) 

CC8=ANTGM*C0S (ZTA (M-1 ) )*DS(M-1 )*US(M-1 ) 

CC9=ANTGN#R(M-l )*DS(M-1 )*US(M-1 ) 

IF ( I 1 1 1 .EQ* 1 ) go TO 83 
IF (M.EQ.MBP) go TO 101 
DNS(M)= (AS(M+1 )-AS(M-1 ) )/(2.*DX) 

GO TO 59 

101 DNS(MBP)=(AS(MBP)-AS(MA2) )/DX 
59 1=M-1 

DNS( I )= (AS( I+l )-AS( I-l ) )/(2.*DX > 
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GO TO 84 

83 DNS<M)=0* 

I =M-1 

DNS( I )=0. 

84 AAI=CC1 
BBl =CC2 

CCI=-(CC8*NS( I )#*2 + CC9*NS( I ) ) 

CCI=CCI+(R( I )+NS( I )*COS{ZTA( !)))♦( (1 .+NS( I )*CK ( I ) )*DS( I )*VS( 1 ) 

1 -DNS ( I )*DS< I )*US( I ) ;*DX 
XX1=CC1 +DX*DS (M )*VS (M )*COS (ZTA (M ) )#CK(M) 

XX2 = CC2 +OX*OS(M )*VS ( M )*R ( M >*CK { M ) 

1 +DX*DS (M )»VS (M )*COS<ZTA (M ) )-COS ( ZTA (M ) >#DS (M> *US 

1 ( M )*DNS (M )»DX 

XX3=-(CC8*(NS (M-l )**2 ) +CC9*NS (M- 1 ) -DX*DS ( M )* VS < M )*R ( M 

1 )+DX*R(M)*DS(M)*US(M)*DNS(M) ) 

NSSl = (-BBl + (BBl**2-4.*AAl*CCI )**0.5 )/ (2**AA1 ) 

NSS2= (-XX2+ (XX2**2-4.*XX1»XX3)**0.5 )/ (2,*XX1 ) 

NSS3= ( NSS2-NSS1 )/NSS2 
331 ADC=0.5 

340 NS<M )=ADC*NSS2+( 1 •-ADC)*NSS1 
IF (NSS3*LT,0, 15) GO TO 82 
75 SS1=NS(M)-NS(M-1 ) 

IF (SSl-GT.O.) GO TO 335 
ADC=ADC+0*0l 1 
GO TO 340 

335 SS2=SS1/NS(M-1 ) 

ZZ=0.02*M-0.02 
IF (SS2.LT.ZZ ) GO TO 82 
ADC=A0C-0«007 
GO TO 340 
82 CONTINUE 
INTGF=0« 

DO 58 N*1 »49 
CC3( 1 »N )=0. 

INTGF=DY (N)/2** { (R(M)+NSS2 *Y(N+1 )*COS(ZTA(M> ) )*D(M*N+1 )#U(M»N+1 >+ 
1 (R(M)+NSS2 *Y(N)*COS(ZTA(M) ) ) *D ( M » N ) *U ( M . N ) ) + INTGF 
CC3(M«N)=INTGF*NSS2 *DS(M) *US ( M ) 

ODD= (CC3<M*N )-CC3 (M-l *N) )/DX 
K = N+1 

EEE= (R( M )+NSS2 *Y ( K ) *COS ( ZTA ( M > ) )*(1 *+NSS2 *Y(K;#CK(M) )*DS(M)*VS(M 
1 )*D(M»K ) 

FFF = - (R (M )+NSS2 *Y (K )*C0S (ZTA (M ) ) )^tDNS (M )*Y (K ) *DS ( M ) *US ( M ) *D ( M* K ) * 
1U(M«K ) 

V ( M « K ) =- ( DDD+FFF ) /EEE 


58 CONTINUE 
ANTGF=0 • 

DO 351 N=1 .49 

ANTGFsDY (N )/2.* ( ( R ( M > +NS ( M ) * Y ( N+ 1 )*COS(ZTA(M) ) )*D(M.N+1 )*U<M.N+1 )+ 
I (R<M)+NS(M)»Y(N)*COS(ZTA(M) ) ) *D ( M . N ) »U ( M . N > )+ANTGF 
CC3(M.N)=ANTGF*NS(M)*DS(M) *US ( M) 

351 continue 
350 P (M.K ) = 1 • 

IF ( I I 1 1 *GT. 1 ) GO TO 60 

DO 62 N 1=1 .50 

N=51-N1 

P(M»N)=P(M.N+1 )-DY<N)*NS(M)*D(M.N)*DS( M ) *U ( M .N ) **2*US < M ) **2*CK ( M 
I ) / (PS ( M )* ( 1 #+NS (M )*Y (N )*CK (M ) ) ) 

62 CONTINUE 
GO TO 28 

60 DVS(M)= (VS(M)-VS(M-l ) )/DX 

DVX(51 > = (V(M.5l )-V(M-l .51 > >/OX 

G8=l •/( 1 *+NS(M)*CK(M ) ) * ( OVS ( M ) /VS (M ) +DVX < 5 1 ) -DNS ( M ) *DVY ( 5 1 )/NS(M) ) 
G9=VS(M )#DVY (51 ) /US ( M ) /NS ( M ) 

G1 0 = -US (M )*CK (M )/(VS(M )♦( 1 • +NS ( M ) *CK ( M ) ) ) 

G7=PS ( M ) / ( DS ( M ) *US ( M ) *VS ( M ) *NS ( M ) ) 

AING1= (G8+G9+G1 0)/G7 
64 DO 68 Nl = 1 .Kl I 
N=51-N1 

IF (M.EQ*MBP) GO TO 93 

DVX(N)= (V(M+1 .N )-V(M-l .N) )/DX/2. 

GO TO 94 

93 DVX(N)= ( V(M0P.N >-V(MA2.N) )/DX 

94 G4=D(M*N)*U(M«N)/ ( 1+NS (M ) *Y ( N 1 *CK ( M ) ) * ( OVS ( M ) ♦ V ( M * N ) /VS ( M ) +OVX ( N ) - 
1DNS(M )*Y (N)*DVY(N )/NS(M) ) 

G5 = VS(M)*D(M»N)*V (M.N)*DVY(N)/US(M >/NS(M) 

G6=-US ( M )*CK (M )*D (M.N )*U(M.N )**2/ ( 1 •+NS(M )*Y (N )*CK (M ) >/VS( M ) 

AInG2= (G4+G5+G6 )/G7 

P(M.N)=P(M.N+1 >+DY(N)* ( A I NGl +A I NG2 ) /2 • 

AING1=A ING2 

68 continue 

28 DO 79 N=1 .Kl 1 
pps=P(M.N)*PS(M) 

PPI=PP*D1NF*VINF**2 

TT=T(M.N)#TS(M) 

TTl =TT»VINF**2/CP I 
105 PP2= (PPl/1013250. )**AL 

IF (TTl •GT.SBOO ) GO TO 71 
CT3=CT»0*9075 



CTT= (CT-CT3 )* (T-4000 )/4500+CT3 
TT2-TTI /CTT 
GO TO 73 
71 TT2=TT1/CT 

73 DD1=0.001292*(PP2/TT2 )*♦< 1 ./AK ) 
61 DD=D01/DINF 

0(M*N )=DO/DS(M ) 

79 CONTINUE 
RETURN 
END 



